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 8104 : FTMachine::FTMachine() : isDryRun(false), image(0), uvwMachine_p(0),
93 8104 : tangentSpecified_p(false), fixMovingSource_p(false),
94 4052 : ephemTableName_p(""),
95 4052 : movingDirShift_p(0.0),
96 4052 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1), romscol_p(nullptr),
97 4052 : useDoubleGrid_p(false),
98 4052 : freqFrameValid_p(false),
99 4052 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
100 4052 : pointingDirCol_p("DIRECTION"),
101 4052 : cfStokes_p(), cfCache_p(), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
102 4052 : canComputeResiduals_p(false), toVis_p(true),
103 20260 : 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 4052 : spectralCoord_p=SpectralCoordinate();
106 4052 : isPseudoI_p=false;
107 4052 : spwChanSelFlag_p=0;
108 4052 : polInUse_p=0;
109 4052 : pop_p = new PolOuterProduct;
110 4052 : ft_p=FFT2D(true);
111 4052 : }
112 :
113 93 : FTMachine::FTMachine(CountedPtr<CFCache>& cfcache,CountedPtr<ConvolutionFunction>& cf):
114 186 : isDryRun(false), image(0), uvwMachine_p(0),
115 186 : tangentSpecified_p(false), fixMovingSource_p(false),
116 93 : ephemTableName_p(""),
117 93 : movingDirShift_p(0.0),
118 93 : distance_p(0.0), lastFieldId_p(-1),lastMSId_p(-1), romscol_p(nullptr),
119 93 : useDoubleGrid_p(false),
120 93 : freqFrameValid_p(false),
121 93 : freqInterpMethod_p(InterpolateArray1D<Double,Complex>::nearestNeighbour),
122 93 : pointingDirCol_p("DIRECTION"),
123 93 : cfStokes_p(), cfCache_p(cfcache), cfs_p(), cfwts_p(), cfs2_p(), cfwts2_p(),
124 93 : convFuncCtor_p(cf),canComputeResiduals_p(false), toVis_p(true), numthreads_p(-1),
125 465 : 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 93 : spectralCoord_p=SpectralCoordinate();
128 93 : isPseudoI_p=false;
129 93 : spwChanSelFlag_p=0;
130 93 : polInUse_p=0;
131 93 : pop_p = new PolOuterProduct;
132 93 : ft_p=FFT2D(true);
133 93 : }
134 :
135 20170 : LogIO& FTMachine::logIO() {return logIO_p;};
136 :
137 : //----------------------------------------------------------------------
138 1132 : FTMachine& FTMachine::operator=(const FTMachine& other)
139 : {
140 1132 : if(this!=&other) {
141 1132 : image=other.image;
142 : //generic selection stuff and state
143 1132 : nAntenna_p=other.nAntenna_p;
144 1132 : distance_p=other.distance_p;
145 1132 : lastFieldId_p=other.lastFieldId_p;
146 1132 : lastMSId_p=other.lastMSId_p;
147 1132 : romscol_p=other.romscol_p;
148 1132 : tangentSpecified_p=other.tangentSpecified_p;
149 1132 : mTangent_p=other.mTangent_p;
150 1132 : mImage_p=other.mImage_p;
151 1132 : mFrame_p=other.mFrame_p;
152 :
153 1132 : nx=other.nx;
154 1132 : ny=other.ny;
155 1132 : npol=other.npol;
156 1132 : nchan=other.nchan;
157 1132 : nvischan=other.nvischan;
158 1132 : nvispol=other.nvispol;
159 1132 : mLocation_p=other.mLocation_p;
160 1132 : if(uvwMachine_p)
161 0 : delete uvwMachine_p;
162 1132 : if(other.uvwMachine_p)
163 0 : uvwMachine_p=new casacore::UVWMachine(*other.uvwMachine_p);
164 : else
165 1132 : uvwMachine_p=0;
166 1132 : doUVWRotation_p=other.doUVWRotation_p;
167 : //Spectral and pol stuff
168 1132 : freqInterpMethod_p=other.freqInterpMethod_p;
169 1132 : spwChanSelFlag_p.resize();
170 1132 : spwChanSelFlag_p=other.spwChanSelFlag_p;
171 1132 : freqFrameValid_p=other.freqFrameValid_p;
172 : //selectedSpw_p.resize();
173 : //selectedSpw_p=other.selectedSpw_p;
174 1132 : imageFreq_p.resize();
175 1132 : imageFreq_p=other.imageFreq_p;
176 1132 : lsrFreq_p.resize();
177 1132 : lsrFreq_p=other.lsrFreq_p;
178 1132 : interpVisFreq_p.resize();
179 1132 : interpVisFreq_p=other.interpVisFreq_p;
180 : //multiChanMap_p=other.multiChanMap_p;
181 1132 : chanMap.resize();
182 1132 : chanMap=other.chanMap;
183 1132 : polMap.resize();
184 1132 : polMap=other.polMap;
185 1132 : nVisChan_p.resize();
186 1132 : nVisChan_p=other.nVisChan_p;
187 1132 : spectralCoord_p=other.spectralCoord_p;
188 1132 : visPolMap_p.resize();
189 1132 : visPolMap_p=other.visPolMap_p;
190 : //doConversion_p.resize();
191 : //doConversion_p=other.doConversion_p;
192 1132 : pointingDirCol_p=other.pointingDirCol_p;
193 : //moving source stuff
194 1132 : movingDir_p=other.movingDir_p;
195 1132 : fixMovingSource_p=other.fixMovingSource_p;
196 1132 : ephemTableName_p = other.ephemTableName_p;
197 1132 : firstMovingDir_p=other.firstMovingDir_p;
198 1132 : movingDirShift_p=other.movingDirShift_p;
199 : //Double precision gridding for those FTMachines that can do
200 1132 : useDoubleGrid_p=other.useDoubleGrid_p;
201 1132 : cfStokes_p = other.cfStokes_p;
202 1132 : polInUse_p = other.polInUse_p;
203 1132 : cfs_p = other.cfs_p;
204 1132 : cfwts_p = other.cfwts_p;
205 1132 : cfs2_p = other.cfs2_p;
206 1132 : cfwts2_p = other.cfwts2_p;
207 1132 : canComputeResiduals_p = other.canComputeResiduals_p;
208 :
209 1132 : pop_p = other.pop_p;
210 1132 : toVis_p = other.toVis_p;
211 1132 : spwFreqSel_p.resize();
212 1132 : spwFreqSel_p = other.spwFreqSel_p;
213 1132 : expandedSpwFreqSel_p = other.expandedSpwFreqSel_p;
214 1132 : expandedSpwConjFreqSel_p = other.expandedSpwConjFreqSel_p;
215 1132 : cmplxImage_p=other.cmplxImage_p;
216 1132 : vbutil_p=other.vbutil_p;
217 1132 : numthreads_p=other.numthreads_p;
218 1132 : pbLimit_p=other.pbLimit_p;
219 1132 : convFuncCtor_p = other.convFuncCtor_p;
220 1132 : sj_p.resize();
221 1132 : sj_p=other.sj_p;
222 1132 : isDryRun=other.isDryRun;
223 1132 : phaseCenterTime_p=other.phaseCenterTime_p;
224 1132 : doneThreadPartition_p=other.doneThreadPartition_p;
225 1132 : xsect_p=other.xsect_p;
226 1132 : ysect_p=other.ysect_p;
227 1132 : nxsect_p=other.nxsect_p;
228 1132 : nysect_p=other.nysect_p;
229 1132 : obsvelconv_p=other.obsvelconv_p;
230 1132 : mtype_p=other.mtype_p;
231 1132 : briggsWeightor_p=other.briggsWeightor_p;
232 1132 : ft_p=other.ft_p;
233 1132 : ftmType_p = other.ftmType_p;
234 1132 : avgPBReady_p = other.avgPBReady_p;
235 : };
236 1132 : 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 3706 : 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 3706 : polInUse_p = 0;
277 3706 : uInt N=0;
278 11630 : for(uInt i=0;i<polMap.nelements();i++) if (polMap(i) > -1) polInUse_p++;
279 3706 : cfStokes_p.resize(polInUse_p);
280 11630 : for(uInt i=0;i<polMap.nelements();i++)
281 7924 : if (polMap(i) > -1) {cfStokes_p(N) = vb.correlationTypes()(i);N++;}
282 : }
283 3706 : }
284 : //----------------------------------------------------------------------
285 3706 : void FTMachine::initMaps(const vi::VisBuffer2& vb) {
286 3706 : logIO() << LogOrigin("FTMachine", "initMaps") << LogIO::NORMAL;
287 :
288 3706 : AlwaysAssert(image, AipsError);
289 :
290 : // Set the frame for the UVWMachine
291 3706 : if(vb.isAttached()){
292 : //mFrame_p=MeasFrame(MEpoch(Quantity(vb.time()(0), "s"), MSColumns(vb.ms()).timeMeas()(0).getRef()), mLocation_p);
293 3706 : if(vbutil_p.null())
294 2299 : vbutil_p=new VisBufferUtil(vb);
295 3706 : romscol_p=new MSColumns(vb.ms());
296 7412 : Unit epochUnit=(romscol_p->time()).keywordSet().asArrayString("QuantumUnits")(IPosition(1,0));
297 3706 : if(!mFrame_p.epoch())
298 1883 : mFrame_p.set(MEpoch(Quantity(vb.time()(0), epochUnit), (romscol_p->timeMeas())(0).getRef()));
299 : else
300 1823 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), epochUnit), (romscol_p->timeMeas())(0).getRef()));
301 3706 : if(!mFrame_p.position())
302 1883 : mFrame_p.set(mLocation_p);
303 : else
304 1823 : mFrame_p.resetPosition(mLocation_p);
305 3706 : if(!mFrame_p.direction())
306 1901 : mFrame_p.set(vbutil_p->getEphemDir(vb, phaseCenterTime_p));
307 : else
308 1805 : mFrame_p.resetDirection(vbutil_p->getEphemDir(vb, phaseCenterTime_p));
309 3706 : }
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 3706 : casacore::CoordinateSystem coords=image->coordinates();
319 3706 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
320 3706 : AlwaysAssert(directionIndex>=0, AipsError);
321 : DirectionCoordinate
322 3706 : directionCoord=coords.directionCoordinate(directionIndex);
323 3706 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
324 3706 : AlwaysAssert(spectralIndex>-1, AipsError);
325 3706 : spectralCoord_p=coords.spectralCoordinate(spectralIndex);
326 :
327 : // get the first position of moving source
328 3706 : 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 30 : mFrame_p.resetEpoch(coords.obsInfo().obsDate());
335 : //Double firstTime=romscol_p->time()(0);
336 :
337 30 : Double firstTime=coords.obsInfo().obsDate().get("s").getValue();
338 : //First convert to HA-DEC or AZEL for parallax correction
339 30 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
340 30 : MDirection tmphadec;
341 30 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
342 : //cerr << "phaseCenterTime_p " << phaseCenterTime_p << endl;
343 10 : tmphadec = MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p > 0.0 ? phaseCenterTime_p : firstTime)), outref1)();
344 20 : MeasComet mcomet(Path((romscol_p->field()).ephemPath(vb.fieldId()(0))).absoluteName());
345 10 : if (mFrame_p.comet())
346 2 : mFrame_p.resetComet(mcomet);
347 : else
348 8 : mFrame_p.set(mcomet);
349 30 : } else if (upcase(movingDir_p.getRefString()).contains("COMET")) {
350 18 : MeasComet mcomet(Path(ephemTableName_p).absoluteName());
351 18 : if (mFrame_p.comet())
352 18 : mFrame_p.resetComet(mcomet);
353 : else
354 0 : mFrame_p.set(mcomet);
355 18 : tmphadec = MDirection::Convert(MDirection(MDirection::COMET), outref1)();
356 18 : } else {
357 2 : tmphadec = MDirection::Convert(movingDir_p, outref1)();
358 : }
359 30 : MDirection::Ref outref(directionCoord.directionType(), mFrame_p);
360 30 : 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 30 : 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 16 : initSourceFreqConv();
376 : }
377 : ///TESTOO
378 : ///waiting for CAS-11060
379 : //firstMovingDir_p=MDirection::Convert(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), outref)();
380 : ////////////////////
381 30 : }
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 3706 : Vector<Double> pixelPhaseCenter(2);
392 3706 : pixelPhaseCenter(0) = Double( image->shape()(0) / 2 );
393 3706 : pixelPhaseCenter(1) = Double( image->shape()(1) / 2 );
394 3706 : directionCoord.toWorld(mImage_p, pixelPhaseCenter);
395 3706 : }
396 :
397 : // Get the object distance in meters
398 3706 : Record info(image->miscInfo());
399 3706 : 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 3706 : initUVWMachine(vb);
409 :
410 3706 : lastFieldId_p=-1;
411 :
412 3706 : lastMSId_p=vb.msId();
413 :
414 : // Set up maps
415 :
416 :
417 :
418 : //Store the image/grid channels freq values
419 : {
420 3706 : Int chanNumbre=image->shape()(3);
421 3706 : Vector<Double> pixindex(chanNumbre);
422 3706 : imageFreq_p.resize(chanNumbre);
423 3706 : Vector<Double> tempStorFreq(chanNumbre);
424 3706 : indgen(pixindex);
425 : // pixindex=pixindex+1.0;
426 19607 : for (Int ll=0; ll< chanNumbre; ++ll){
427 15901 : if( !spectralCoord_p.toWorld(tempStorFreq(ll), pixindex(ll))){
428 0 : logIO() << "Cannot get imageFreq " << LogIO::EXCEPTION;
429 :
430 : }
431 : }
432 3706 : convertArray(imageFreq_p,tempStorFreq);
433 3706 : }
434 : //Destroy any conversion layer Freq coord if freqframe is not valid
435 3706 : if(!freqFrameValid_p){
436 846 : MFrequency::Types imageFreqType=spectralCoord_p.frequencySystem();
437 846 : spectralCoord_p.setFrequencySystem(imageFreqType);
438 846 : spectralCoord_p.setReferenceConversion(imageFreqType,
439 1692 : MEpoch(Quantity(vb.time()(0), "s")),
440 846 : mLocation_p,
441 846 : 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 3706 : nvischan = vb.getFrequencies(0).nelements();
451 3706 : interpVisFreq_p.resize();
452 3706 : interpVisFreq_p=vb.getFrequencies(0);
453 :
454 : // Polarization map
455 3706 : visPolMap_p.resize();
456 3706 : 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 3706 : chanMap.resize();
461 3706 : matchChannel(vb);
462 : //chanMap=multiChanMap_p[vb.spectralWindows()(0)];
463 3706 : 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 3706 : if(max(chanMap)>=nchan||min(chanMap)<-2) {
471 0 : logIO() << "Illegal Channel Map: " << chanMap << LogIO::EXCEPTION;
472 : }
473 :
474 :
475 3706 : initPolInfo(vb);
476 3706 : Vector<Int> intpolmap(visPolMap_p.nelements());
477 11630 : for (uInt kk=0; kk < intpolmap.nelements(); ++kk){
478 7924 : intpolmap[kk]=Int(visPolMap_p[kk]);
479 : }
480 3706 : pop_p->initCFMaps(intpolmap, polMap);
481 :
482 :
483 :
484 :
485 :
486 :
487 :
488 3706 : }
489 :
490 3569 : 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 3569 : casacore::CoordinateSystem coords=image->coordinates();
495 3569 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
496 3569 : AlwaysAssert(directionIndex>=0, AipsError);
497 3569 : auto const directionCoord=coords.directionCoordinate(directionIndex);
498 7138 : Vector<Double> equal= (mImage_p.getAngle()-
499 10707 : vbutil_p->getPhaseCenter(vb, phaseCenterTime_p).getAngle()).getValue();
500 10707 : if((abs(equal(0)) < abs(directionCoord.increment()(0)))
501 10707 : && (abs(equal(1)) < abs(directionCoord.increment()(1)))){
502 3020 : doUVWRotation_p=false;
503 : } else {
504 549 : doUVWRotation_p=true;
505 : }
506 :
507 3569 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
508 3569 : String observatory;
509 3569 : if(vb.isAttached())
510 3569 : observatory=(vb.subtableColumns().observation()).telescopeName()(0);
511 : else
512 0 : throw(AipsError("Cannot define frame because of no access to OBSERVATION table"));
513 7138 : if(observatory.contains("ATCA") || observatory.contains("DRAO")
514 7138 : || 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 7138 : uvwMachine_p=new casacore::UVWMachine(mImage_p, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
519 3569 : false, tangentSpecified_p);
520 : }
521 3569 : AlwaysAssert(uvwMachine_p, AipsError);
522 :
523 3569 : phaseShifter_p=new UVWMachine(*uvwMachine_p);
524 3569 : }
525 :
526 2773 : void FTMachine::initBriggsWeightor(vi::VisibilityIterator2& vi){
527 : ///Lastly initialized Briggs cube weighting scheme
528 2773 : if(!briggsWeightor_p.null()){
529 29 : String error;
530 29 : Record rec;
531 29 : AlwaysAssert(image, AipsError);
532 29 : if(!toRecord(error, rec))
533 0 : throw (AipsError("Could not initialize BriggsWeightor"));
534 29 : String wgtcolname=briggsWeightor_p->initImgWeightCol(vi, *image, rec);
535 29 : tempFileNames_p.resize(tempFileNames_p.nelements()+1, True);
536 29 : tempFileNames_p[tempFileNames_p.nelements()-1]=wgtcolname;
537 :
538 29 : }
539 2773 : }
540 :
541 4142 : FTMachine::~FTMachine()
542 : {
543 4142 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
544 4142 : }
545 :
546 :
547 16 : void FTMachine::initSourceFreqConv(){
548 16 : MRadialVelocity::Types refvel=MRadialVelocity::GEO;
549 16 : if(mFrame_p.comet()){
550 : //Has a ephem table
551 16 : if(((mFrame_p.comet())->getTopo().getLength("km").getValue()) > 1.0e-3){
552 11 : 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 32 : obsvelconv_p=MRadialVelocity::Convert (MRadialVelocity(MVRadialVelocity(0.0),
602 32 : MRadialVelocity::Ref(MRadialVelocity::TOPO, mFrame_p)),
603 48 : MRadialVelocity::Ref(refvel));
604 :
605 16 : }
606 :
607 :
608 2768 : Long FTMachine::estimateRAM(const CountedPtr<SIImageStore>& imstor){
609 : //not set up yet
610 2768 : if(!image && !imstor)
611 0 : return -1;
612 2768 : Long npixels=0;
613 2768 : if(image)
614 14 : npixels=((image->shape()).product())/1024;
615 : else{
616 2754 : if((imstor->getShape()).product() !=0)
617 2614 : npixels=(imstor->getShape()).product()/1024;
618 : }
619 2768 : if(npixels==0) npixels=1; //1 kPixels is minimum then
620 2768 : Long factor=sizeof(Complex);
621 2768 : if(!toVis_p && useDoubleGrid_p)
622 0 : factor=sizeof(DComplex);
623 2768 : return (npixels*factor);
624 : }
625 :
626 2393 : void FTMachine::shiftFreqToSource(Vector<Double>& freqs){
627 2393 : MDoppler dopshift;
628 2393 : MEpoch ep(mFrame_p.epoch());
629 2393 : 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 2393 : MEpoch::Convert toUT(ep, MEpoch::UT);
633 2393 : MVRadialVelocity cometvel;
634 2393 : (*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 2393 : dopshift=MDoppler(Quantity(-cometvel.get("km/s").getValue("km/s")+obsvelconv_p().get("km/s").getValue("km/s") , "km/s"), MDoppler::RELATIVISTIC);
639 :
640 2393 : }
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 2393 : Vector<Double> newfreqs=dopshift.shiftFrequency(freqs);
660 2393 : freqs=newfreqs;
661 2393 : }
662 :
663 770555 : 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 770555 : Cube<Complex> origdata;
670 770555 : Cube<Bool> modflagCube;
671 : // Read flags from the vb.
672 770555 : setSpectralFlag(vb,modflagCube);
673 770555 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
674 : //if(doConversion_p[vb.spectralWindows()[0]]){
675 770555 : if(freqFrameValid_p){
676 602953 : visFreq.resize(lsrFreq_p.shape());
677 602953 : convertArray(visFreq, lsrFreq_p);
678 : }
679 : else{
680 167602 : convertArray(visFreq, vb.getFrequencies(0));
681 167602 : lsrFreq_p.resize();
682 167602 : lsrFreq_p=vb.getFrequencies(0);
683 : }
684 770555 : if(type==FTMachine::MODEL){
685 0 : origdata.reference(vb.visCubeModel());
686 : }
687 770555 : else if(type==FTMachine::CORRECTED){
688 269453 : origdata.reference(vb.visCubeCorrected());
689 : }
690 501102 : else if(type==FTMachine::OBSERVED){
691 276023 : origdata.reference(vb.visCube());
692 : }
693 225079 : else if(type==FTMachine::PSF){
694 : // make sure its a size 0 data ...psf
695 : //so avoid reading any data from disk
696 225079 : origdata.resize();
697 :
698 : }
699 : else{
700 0 : throw(AipsError("Don't know which column is being regridded"));
701 : }
702 770555 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour) || (vb.nChannels()==1)){
703 568627 : data.reference(origdata);
704 : // do something here for apply flag based on spw chan sels
705 : // e.g.
706 :
707 :
708 568627 : flags.resize(modflagCube.shape());
709 568627 : flags=0;
710 : //flags(vb.flagCube())=true;
711 :
712 568627 : flags(modflagCube)=true;
713 :
714 568627 : weight.reference(wt);
715 568627 : interpVisFreq_p.resize();
716 568627 : interpVisFreq_p=lsrFreq_p;
717 : //cerr << "INTERPTOGRID " << interpVisFreq_p.nelements() << " vb.nchan " << vb.nChannels() << endl;
718 568627 : return false;
719 : }
720 :
721 201928 : Cube<Bool>flag;
722 :
723 : //okay at this stage we have at least 2 channels
724 201928 : 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 201928 : if(width > 1.0){
734 50846 : Double minVF=min(visFreq);
735 50846 : Double maxVF=max(visFreq);
736 50846 : Double minIF=min(imageFreq_p);
737 50846 : Double maxIF=max(imageFreq_p);
738 101692 : if( ((minIF-fabs(imageFreq_p[1]-imageFreq_p[0])/2.0) > maxVF) ||
739 50846 : ((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 50846 : uInt where=0;
750 : //Double interpwidth=visFreq[1]-visFreq[0];
751 50846 : 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 50846 : if(minIF < minVF){ // Need to find the first image-channel with data in it
755 5946 : where=binarySearchBrackets(found, imageFreq_p, minVF, imageFreq_p.nelements());
756 5946 : if(where != imageFreq_p.nelements()){
757 5946 : minIF=imageFreq_p[where];
758 : }
759 : }
760 :
761 50846 : if(maxIF > maxVF){
762 4409 : where=binarySearchBrackets(found, imageFreq_p, maxVF, imageFreq_p.nelements());
763 4409 : if(where!= imageFreq_p.nelements()){
764 4409 : 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 50846 : Int ninterpchan=(Int)ceil((maxIF-minIF+fabs(imageFreq_p[1]-imageFreq_p[0]))/fabs(interpwidth))+2;
772 50846 : chanMap.resize(ninterpchan);
773 50846 : chanMap.set(-1);
774 50846 : interpVisFreq_p.resize(ninterpchan);
775 50846 : interpVisFreq_p[0]=(interpwidth > 0) ? minIF : maxIF;
776 50846 : if(freqInterpMethod_p==InterpolateArray1D<Double, Complex>::linear)
777 50846 : interpVisFreq_p[0]-=interpwidth;
778 50846 : if(freqInterpMethod_p==InterpolateArray1D<Double, Complex>::cubic)
779 0 : interpVisFreq_p[0]-=2.0*interpwidth;
780 50846 : Double startedge=abs(imageFreq_p[1]-imageFreq_p[0])/2.0 -abs(interpwidth)/2.0;
781 50846 : interpVisFreq_p[0] =(interpwidth >0) ? (interpVisFreq_p[0]-startedge):(interpVisFreq_p[0]+startedge);
782 :
783 12586638 : for (Int k=1; k < ninterpchan; ++k){
784 12535792 : interpVisFreq_p[k] = interpVisFreq_p[k-1]+ interpwidth;
785 : }
786 50846 : Double halfdiff=fabs((imageFreq_p[1]-imageFreq_p[0])/2.0);
787 12637484 : 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 12586638 : Int which=-1;
792 509498740 : 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 496912102 : if( (interpVisFreq_p[k] >= (imageFreq_p[j]-halfdiff)) && (interpVisFreq_p[k] < (imageFreq_p[j]+halfdiff)))
795 12489335 : which=j;
796 : }
797 12586638 : if((which > -1) && (which < Int(imageFreq_p.nelements()))){
798 12486935 : 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 151082 : interpVisFreq_p.resize(imageFreq_p.nelements());
815 151082 : convertArray(interpVisFreq_p, imageFreq_p);
816 151082 : chanMap.resize(interpVisFreq_p.nelements());
817 151082 : indgen(chanMap);
818 : }
819 201928 : 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 143186 : Cube<Complex> flipdata;
823 143186 : 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 143186 : swapyz(flipflag,modflagCube);
830 143186 : swapyz(flipdata,origdata);
831 : InterpolateArray1D<Double,Complex>::
832 143186 : interpolate(data,flag,interpVisFreq_p,visFreq,flipdata,flipflag,freqInterpMethod_p, False, False);
833 143186 : flipdata.resize();
834 143186 : swapyz(flipdata,data);
835 143186 : data.resize();
836 143186 : data.reference(flipdata);
837 143186 : flipflag.resize();
838 143186 : swapyz(flipflag,flag);
839 143186 : flag.resize();
840 143186 : flag.reference(flipflag);
841 : // Note : 'flag' will get augmented with the flags coming out of weight interpolation
842 143186 : }
843 : else
844 : { // get the flag array to the correct shape.
845 : // This will get filled at the end of weight-interpolation.
846 58742 : flag.resize(vb.nCorrelations(), interpVisFreq_p.nelements(), vb.nRows());
847 58742 : 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 201928 : InterpolateArray1D<casacore::Double,casacore::Complex>::InterpolationMethod weightinterp=freqInterpMethod_p;
854 :
855 201928 : if(!briggsWeightor_p.null()){
856 9190 : weightinterp= InterpolateArray1D<casacore::Double,casacore::Complex>::nearestNeighbour;
857 : }
858 : //InterpolateArray1D<casacore::Double,casacore::Complex>::InterpolationMethod weightinterp=InterpolateArray1D<casacore::Double,casacore::Complex>::nearestNeighbour;
859 201928 : Matrix<Bool> chanflag(wt.shape());
860 201928 : AlwaysAssert( chanflag.shape()[0]==modflagCube.shape()[1], AipsError);
861 201928 : AlwaysAssert( chanflag.shape()[1]==modflagCube.shape()[2], AipsError);
862 201928 : chanflag=false;
863 605916 : for(uInt pol=0;pol<modflagCube.shape()[0];pol++)
864 403988 : 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 201928 : Matrix<Float> flipweight;
871 201928 : flipweight=transpose(wt);
872 201928 : Matrix<Bool> flipchanflag;
873 201928 : flipchanflag=transpose(chanflag);
874 201928 : Matrix<Bool> tempoutputflag;
875 : InterpolateArray1D<Double,Float>::
876 201928 : interpolate(weight,tempoutputflag, interpVisFreq_p, visFreq,flipweight,flipchanflag,weightinterp, False, False);
877 201928 : flipweight.resize();
878 201928 : flipweight=transpose(weight);
879 201928 : weight.resize();
880 201928 : weight.reference(flipweight);
881 201928 : flipchanflag.resize();
882 201928 : flipchanflag=transpose(tempoutputflag);
883 201928 : tempoutputflag.resize();
884 201928 : 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 201928 : AlwaysAssert( tempoutputflag.shape()[0]==flag.shape()[1], AipsError);
896 201928 : AlwaysAssert( tempoutputflag.shape()[1]==flag.shape()[2], AipsError);
897 605916 : for(uInt pol=0;pol<flag.shape()[0];pol++)
898 403988 : flag.yzPlane(pol) = tempoutputflag | flag.yzPlane(pol);
899 :
900 : // Fill the output array of image-channel flags.
901 201928 : flags.resize(flag.shape());
902 201928 : flags=0;
903 201928 : flags(flag)=true;
904 :
905 201928 : return true;
906 770555 : }
907 :
908 187791 : void FTMachine::getInterpolateArrays(const vi::VisBuffer2& vb,
909 : Cube<Complex>& data, Cube<Int>& flags){
910 :
911 187791 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
912 :
913 : //if(doConversion_p[vb.spectralWindows()[0]]){
914 187791 : if(freqFrameValid_p){
915 140289 : convertArray(visFreq, lsrFreq_p);
916 : }
917 : else{
918 47502 : convertArray(visFreq, vb.getFrequencies(0));
919 : }
920 :
921 :
922 :
923 187791 : if((imageFreq_p.nelements()==1) || (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour)|| (vb.nChannels()==1) ){
924 134197 : Cube<Bool> modflagCube;
925 134197 : setSpectralFlag(vb,modflagCube);
926 :
927 134197 : data.reference(vb.visCubeModel());
928 : //flags.resize(vb.flagCube().shape());
929 134197 : flags.resize(modflagCube.shape());
930 134197 : flags=0;
931 : //flags(vb.flagCube())=true;
932 134197 : flags(modflagCube)=true;
933 134197 : interpVisFreq_p.resize();
934 134197 : interpVisFreq_p=visFreq;
935 134197 : return;
936 134197 : }
937 :
938 53594 : data.resize(vb.nCorrelations(), imageFreq_p.nelements(), vb.nRows());
939 53594 : flags.resize(vb.nCorrelations(), imageFreq_p.nelements(), vb.nRows());
940 53594 : data.set(Complex(0.0,0.0));
941 53594 : flags.set(0);
942 : //no need to degrid channels that does map over this vb
943 53594 : Int maxchan=max(chanMap);
944 1008007 : for (uInt k =0 ; k < chanMap.nelements() ; ++k){
945 954413 : if(chanMap(k)==-1)
946 147324 : chanMap(k)=maxchan;
947 : }
948 53594 : Int minchan=min(chanMap);
949 53594 : if(minchan==maxchan)
950 0 : minchan=-1;
951 :
952 :
953 53594 : for(Int k = 0; k < minchan; ++k)
954 0 : flags.xzPlane(k).set(1);
955 :
956 73274 : for(uInt k = maxchan + 1; k < imageFreq_p.nelements(); ++k)
957 19680 : flags.xzPlane(k).set(1);
958 :
959 53594 : interpVisFreq_p.resize(imageFreq_p.nelements());
960 53594 : convertArray(interpVisFreq_p, imageFreq_p);
961 53594 : chanMap.resize(imageFreq_p.nelements());
962 53594 : indgen(chanMap);
963 187791 : }
964 :
965 187791 : Bool FTMachine::interpolateFrequencyFromgrid(vi::VisBuffer2& vb,
966 : Cube<Complex>& data,
967 : FTMachine::Type type){
968 :
969 : Cube<Complex> *origdata;
970 187791 : Vector<Double> visFreq(vb.getFrequencies(0).nelements());
971 :
972 : //if(doConversion_p[vb.spectralWindows()[0]]){
973 187791 : if(freqFrameValid_p){
974 140289 : convertArray(visFreq, lsrFreq_p);
975 : }
976 : else{
977 47502 : convertArray(visFreq, vb.getFrequencies(0));
978 : }
979 :
980 187791 : if(type==FTMachine::MODEL){
981 187791 : 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 187791 : Double width=fabs(imageFreq_p[1]-imageFreq_p[0])/fabs(visFreq[1]-visFreq[0]);
996 :
997 187791 : if((imageFreq_p.nelements()==1) ||
998 255777 : (vb.nChannels()==1) ||
999 67986 : (freqInterpMethod_p== InterpolateArray1D<Double, Complex>::nearestNeighbour) ){
1000 134197 : interpVisFreq_p=visFreq;
1001 : //cerr << "INTERPFROMGRID " << interpVisFreq_p << " vb.nchan " << vb.nChannels() << endl;
1002 134197 : origdata->reference(data);
1003 134197 : interpVisFreq_p=visFreq;
1004 134197 : 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 53594 : Cube<Complex> flipgrid;
1010 53594 : flipgrid.resize();
1011 53594 : swapyz(flipgrid,data);
1012 53594 : Vector<Double> newImFreq;
1013 53594 : 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 53594 : if(width > 1.0){
1019 4946 : Int newNchan=Int(std::floor(width))*imageFreq_p.nelements();
1020 4946 : newImFreq.resize(newNchan);
1021 4946 : Double newIncr= (imageFreq_p[1]-imageFreq_p[0])/std::floor(width);
1022 4946 : Double newStart=imageFreq_p[0]-(imageFreq_p[1]-imageFreq_p[0])/2.0+newIncr/2.0;
1023 4946 : Cube<Complex> newflipgrid(flipgrid.shape()[0], flipgrid.shape()[1], newNchan);
1024 :
1025 59408 : for (Int k=0; k < newNchan; ++k){
1026 54462 : newImFreq[k]=newStart+k*newIncr;
1027 54462 : Int oldchan=k/Int(std::floor(width));
1028 54462 : 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 4946 : flipgrid.resize();
1036 4946 : flipgrid.reference(newflipgrid);
1037 :
1038 4946 : }
1039 107188 : Cube<Complex> flipdata((origdata->shape())(0),(origdata->shape())(2),
1040 160782 : (origdata->shape())(1)) ;
1041 53594 : 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 53594 : interpolate(flipdata,visFreq, newImFreq, flipgrid,freqInterpMethod_p);
1054 :
1055 :
1056 :
1057 53594 : Cube<Bool> copyOfFlag;
1058 : //Vector<Int> mychanmap=multiChanMap_p[vb.spectralWindows()[0]];
1059 53594 : matchChannel(vb);
1060 53594 : copyOfFlag.assign(vb.flagCube());
1061 1008007 : for (uInt k=0; k< chanMap.nelements(); ++ k)
1062 954413 : if(chanMap(k) == -1)
1063 147324 : copyOfFlag.xzPlane(k).set(true);
1064 53594 : flipgrid.resize();
1065 53594 : swapyz(flipgrid, copyOfFlag, flipdata);
1066 : //swapyz(flipgrid,flipdata);
1067 53594 : vb.setVisCubeModel(flipgrid);
1068 :
1069 53594 : return true;
1070 187791 : }
1071 :
1072 :
1073 112029 : 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 112029 : if(lastMSId_p != vb.msId())
1083 2 : romscol_p=new MSColumns(vb.ms());
1084 112029 : if((vb.fieldId()(0)!=lastFieldId_p) || (vb.msId()!=lastMSId_p)){
1085 2069 : doUVWRotation_p=true;
1086 : }
1087 : else{
1088 : //if above failed it still can be changing if polynome phasecenter or ephem
1089 :
1090 109960 : if( (vb.subtableColumns().field().numPoly()(lastFieldId_p) >0) || (! (vb.subtableColumns().field().ephemerisId().isNull()) && (vb.subtableColumns().field().ephemerisId()(lastFieldId_p) > -1)))
1091 1064 : doUVWRotation_p=True;
1092 : }
1093 112029 : if(doUVWRotation_p || fixMovingSource_p){
1094 :
1095 112029 : mFrame_p.epoch() != 0 ?
1096 224058 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
1097 112029 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), (romscol_p->timeMeas())(0).getRef()));
1098 : MDirection::Types outType;
1099 112029 : MDirection::getType(outType, mImage_p.getRefString());
1100 224058 : MDirection phasecenter=MDirection::Convert(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), MDirection::Ref(outType, mFrame_p))();
1101 112029 : MDirection inFieldPhaseCenter=phasecenter;
1102 :
1103 112029 : if(fixMovingSource_p){
1104 :
1105 : //First convert to HA-DEC or AZEL for parallax correction
1106 1120 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1107 1120 : MDirection tmphadec;
1108 1120 : if(upcase(movingDir_p.getRefString()).contains("APP")){
1109 560 : tmphadec=MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p)), outref1)();
1110 : }
1111 : else{
1112 560 : tmphadec=MDirection::Convert(movingDir_p, outref1)();
1113 : }
1114 1120 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
1115 1120 : 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 1120 : movingDirShift_p=MVDirection(sourcenow.getAngle()-firstMovingDir_p.getAngle());
1119 : // cerr << "shift " << movingDirShift_p.getAngle() << endl;
1120 1120 : inFieldPhaseCenter.shift(movingDirShift_p, False);
1121 1120 : }
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 112029 : if(doUVWRotation_p || fixMovingSource_p) {
1128 :
1129 112029 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
1130 112029 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
1131 112029 : 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 224058 : uvwMachine_p=new UVWMachine(inFieldPhaseCenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1141 112029 : false, false);
1142 112029 : phaseShifter_p=new UVWMachine(mImage_p, phasecenter, mFrame_p,
1143 112029 : false, false);
1144 : }
1145 112029 : }
1146 :
1147 112029 : lastFieldId_p=vb.fieldId()(0);
1148 112029 : lastMSId_p=vb.msId();
1149 :
1150 :
1151 112029 : AlwaysAssert(uvwMachine_p, AipsError);
1152 :
1153 : // Always force a recalculation
1154 112029 : uvwMachine_p->reCalculate();
1155 112029 : phaseShifter_p->reCalculate();
1156 :
1157 : // Now do the conversions
1158 112029 : uInt nrows=dphase.nelements();
1159 112029 : Vector<Double> thisRow(3);
1160 112029 : 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 21150174 : for (uInt irow=0; irow<nrows;++irow) {
1170 21038145 : thisRow.assign(uvw.column(irow));
1171 : //cerr << " uvw " << thisRow ;
1172 : // This is for frame change
1173 21038145 : uvwMachine_p->convertUVW(dphase(irow), thisRow);
1174 : // This is for correlator phase center change
1175 21038145 : 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 21038145 : RotMatrix rotMat=phaseShifter_p->rotationUVW();
1184 21038145 : uvw.column(irow)(0)=thisRow(0)*rotMat(0,0)+thisRow(1)*rotMat(1,0);
1185 21038145 : uvw.column(irow)(1)=thisRow(1)*rotMat(1,1)+thisRow(0)*rotMat(0,1);
1186 21038145 : uvw.column(irow)(2)=thisRow(0)*rotMat(0,2)+thisRow(1)*rotMat(1,2)+thisRow(2)*rotMat(2,2);
1187 21038145 : dphase(irow)+= rotphase(0)*uvw.column(irow)(0)+rotphase(1)*uvw.column(irow)(1);
1188 21038145 : }
1189 :
1190 :
1191 112029 : }
1192 112029 : }
1193 :
1194 :
1195 677091 : void FTMachine::rotateUVW(Matrix<Double>& uvw, Vector<Double>& dphase,
1196 : const vi::VisBuffer2& vb)
1197 : {
1198 :
1199 677091 : if(lastMSId_p != vb.msId())
1200 86 : 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 677091 : if((vb.fieldId()(0)!=lastFieldId_p) || (vb.msId()!=lastMSId_p)){
1206 3119 : doUVWRotation_p=true;
1207 :
1208 : }
1209 : else{
1210 : //if above failed it still can be changing if polynome phasecenter or ephem
1211 673972 : if( (vb.subtableColumns().field().numPoly()(lastFieldId_p) >0) || (! (vb.subtableColumns().field().ephemerisId().isNull()) &&(vb.subtableColumns().field().ephemerisId()(lastFieldId_p) > -1)))
1212 323 : doUVWRotation_p=True;
1213 :
1214 : }
1215 677091 : if(doUVWRotation_p || tangentSpecified_p || fixMovingSource_p){
1216 677091 : ok();
1217 :
1218 677091 : mFrame_p.epoch() != 0 ?
1219 1354182 : mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(0), "s"))):
1220 :
1221 677091 : mFrame_p.set(mLocation_p, MEpoch(Quantity(vb.time()(0), "s"), (romscol_p->timeMeas())(0).getRef()));
1222 :
1223 677091 : MDirection phasecenter=mImage_p;
1224 677091 : if(fixMovingSource_p){
1225 :
1226 : //First convert to HA-DEC or AZEL for parallax correction
1227 340 : MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1228 340 : MDirection tmphadec;
1229 340 : if(upcase(movingDir_p.getRefString()).contains("APP")){
1230 80 : tmphadec=MDirection::Convert((vbutil_p->getEphemDir(vb, phaseCenterTime_p)), outref1)();
1231 : }
1232 : else{
1233 260 : 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 340 : MDirection::Ref outref(mImage_p.getRef().getType(), mFrame_p);
1239 340 : 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 340 : movingDirShift_p=MVDirection(sourcenow.getAngle()-firstMovingDir_p.getAngle());
1243 340 : phasecenter.shift(movingDirShift_p, False);
1244 : //cerr << sourcenow.toString() <<" "<<(vbutil_p->getPhaseCenter(vb, phaseCenterTime_p)).toString() << " difference " << firstMovingDir_p.getAngle() - sourcenow.getAngle() << endl;
1245 340 : }
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 677091 : if(doUVWRotation_p || fixMovingSource_p) {
1252 :
1253 677091 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
1254 677091 : if(uvwMachine_p) delete uvwMachine_p; uvwMachine_p=0;
1255 677091 : 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 1354182 : uvwMachine_p=new UVWMachine(phasecenter, vbutil_p->getPhaseCenter(vb, phaseCenterTime_p), mFrame_p,
1263 677091 : false,tangentSpecified_p);
1264 : }
1265 677091 : }
1266 :
1267 677091 : lastFieldId_p=vb.fieldId()(0);
1268 677091 : lastMSId_p=vb.msId();
1269 :
1270 :
1271 677091 : AlwaysAssert(uvwMachine_p, AipsError);
1272 :
1273 : // Always force a recalculation
1274 677091 : uvwMachine_p->reCalculate();
1275 :
1276 : // Now do the conversions
1277 677091 : uInt nrows=dphase.nelements();
1278 677091 : Vector<Double> thisRow(3);
1279 677091 : thisRow=0.0;
1280 : uInt irow;
1281 : //#pragma omp parallel default(shared) private(irow,thisRow)
1282 : {
1283 : //#pragma omp for
1284 231744562 : for (irow=0; irow<nrows;++irow) {
1285 231067471 : thisRow.reference(uvw.column(irow));
1286 231067471 : convUVW(dphase(irow), thisRow);
1287 : }
1288 :
1289 : }//end pragma
1290 677091 : }
1291 677091 : }
1292 231067471 : void FTMachine::convUVW(Double& dphase, Vector<Double>& thisrow){
1293 : //for (uInt i=0;i<3;i++) thisRow(i)=uvw(i,row);
1294 231067471 : uvwMachine_p->convertUVW(dphase, thisrow);
1295 : //for (uint i=0;i<3;i++) uvw(i,row)=thisRow(i);
1296 :
1297 231067471 : }
1298 :
1299 :
1300 70200 : 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 70200 : Int rowoff=row*nvchan;
1305 : Double phase;
1306 : Double pos;
1307 70200 : Int nel= doW ? 3 : 2;
1308 7090200 : for (Int f=0; f<nvchan; ++f){
1309 21060000 : for (Int k=0; k <2; ++k){
1310 14040000 : pos=(scale[k])*uvw[3*row+k]*(freq[f])/C::c+((offset[k])+1.0);
1311 14040000 : loc[(rowoff+f)*nel+k]=std::lround(pos);
1312 14040000 : 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 7020000 : phase=-Double(2.0)*C::pi*dphase[row]*(freq[f])/C::c;
1316 7020000 : phasor[rowoff+f]=Complex(cos(phase), sin(phase));
1317 :
1318 : ///This is for W-Projection
1319 7020000 : 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 70200 : }
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 789120 : 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 789120 : ok();
1349 :
1350 789120 : 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 789120 : }
1416 :
1417 0 : void FTMachine::ok() {
1418 0 : AlwaysAssert(image, AipsError);
1419 0 : AlwaysAssert(uvwMachine_p, AipsError);
1420 0 : }
1421 :
1422 48 : Bool FTMachine::toRecord(String& error, RecordInterface& outRecord,
1423 : Bool withImage, const String diskimage) {
1424 : // Save the FTMachine to a Record
1425 : //
1426 48 : outRecord.define("name", this->name());
1427 48 : if(withImage){
1428 19 : if(image==nullptr)
1429 0 : throw(AipsError("Programmer error: saving to record without proper initialization"));
1430 19 : CoordinateSystem cs=image->coordinates();
1431 19 : DirectionCoordinate dircoord=cs.directionCoordinate(cs.findCoordinate(Coordinate::DIRECTION));
1432 19 : dircoord.setReferenceValue(mImage_p.getAngle().getValue());
1433 19 : if(diskimage != ""){
1434 : try{
1435 38 : PagedImage<Complex> imCopy(TiledShape(toVis_p ? griddedData.shape(): image->shape()), image->coordinates(), diskimage);
1436 19 : toVis_p ? imCopy.put(griddedData) : imCopy.copyData(*image);
1437 19 : ImageUtilities::copyMiscellaneous(imCopy, *image);
1438 19 : Vector<Double> pixcen(2);
1439 19 : pixcen(0)=Double(imCopy.shape()(0)/2); pixcen(1)=Double(imCopy.shape()(1)/2);
1440 19 : dircoord.setReferencePixel(pixcen);
1441 19 : cs.replaceCoordinate(dircoord, cs.findCoordinate(Coordinate::DIRECTION));
1442 19 : imCopy.setCoordinateInfo(cs);
1443 19 : }
1444 0 : catch(...){
1445 0 : throw(AipsError(String("Failed to save model image "+diskimage+String(" to disk"))));
1446 0 : }
1447 19 : 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 19 : }
1463 48 : outRecord.define("nantenna", nAntenna_p);
1464 48 : outRecord.define("distance", distance_p);
1465 48 : outRecord.define("lastfieldid", lastFieldId_p);
1466 48 : outRecord.define("lastmsid", lastMSId_p);
1467 48 : outRecord.define("tangentspecified", tangentSpecified_p);
1468 48 : saveMeasure(outRecord, String("mtangent_rec"), error, mTangent_p);
1469 48 : saveMeasure(outRecord, "mimage_rec", error, mImage_p);
1470 : //mFrame_p not necessary to save as it is generated from mLocation_p
1471 48 : outRecord.define("nx", nx);
1472 48 : outRecord.define("ny", ny);
1473 48 : outRecord.define("npol", npol);
1474 48 : outRecord.define("nchan", nchan);
1475 48 : outRecord.define("nvischan", nvischan);
1476 48 : outRecord.define("nvispol", nvispol);
1477 : //no need to save uvwMachine_p
1478 48 : outRecord.define("douvwrotation", doUVWRotation_p);
1479 48 : outRecord.define("freqinterpmethod", static_cast<Int>(freqInterpMethod_p));
1480 48 : outRecord.define("spwchanselflag", spwChanSelFlag_p);
1481 48 : outRecord.define("freqframevalid", freqFrameValid_p);
1482 : //outRecord.define("selectedspw", selectedSpw_p);
1483 48 : outRecord.define("imagefreq", imageFreq_p);
1484 48 : outRecord.define("lsrfreq", lsrFreq_p);
1485 48 : outRecord.define("interpvisfreq", interpVisFreq_p);
1486 48 : 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 48 : outRecord.define("chanmap", chanMap);
1491 48 : outRecord.define("polmap", polMap);
1492 48 : outRecord.define("nvischanmulti", nVisChan_p);
1493 : //save moving source related variables
1494 48 : storeMovingSourceState(error, outRecord);
1495 : //outRecord.define("doconversion", doConversion_p);
1496 48 : outRecord.define("pointingdircol", pointingDirCol_p);
1497 48 : outRecord.define("usedoublegrid", useDoubleGrid_p);
1498 48 : outRecord.define("cfstokes", cfStokes_p);
1499 48 : outRecord.define("polinuse", polInUse_p);
1500 48 : outRecord.define("tovis", toVis_p);
1501 48 : outRecord.define("sumweight", sumWeight);
1502 48 : outRecord.define("numthreads", numthreads_p);
1503 48 : outRecord.define("phasecentertime", phaseCenterTime_p);
1504 : //Need to serialized sj_p...the user has to set the sj_p after recovering from record
1505 48 : return true;
1506 48 : };
1507 :
1508 240 : Bool FTMachine::saveMeasure(RecordInterface& rec, const String& name, String& err, const Measure& meas){
1509 240 : Record tmprec;
1510 240 : MeasureHolder mh(meas);
1511 240 : if(mh.toRecord(err, tmprec)){
1512 240 : rec.defineRecord(name, tmprec);
1513 240 : return true;
1514 : }
1515 0 : return false;
1516 240 : }
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 48 : Bool FTMachine::storeMovingSourceState(String& error, RecordInterface& outRecord){
1651 :
1652 48 : Bool retval=True;
1653 48 : retval=retval && saveMeasure(outRecord, "mlocation_rec", error, mLocation_p);
1654 48 : spectralCoord_p.save(outRecord, "spectralcoord");
1655 48 : retval=retval && saveMeasure(outRecord, "movingdir_rec", error, movingDir_p);
1656 48 : outRecord.define("fixmovingsource", fixMovingSource_p);
1657 48 : retval=retval && saveMeasure(outRecord, "firstmovingdir_rec", error, firstMovingDir_p);
1658 48 : movingDirShift_p=MVDirection(0.0);
1659 48 : if( mFrame_p.comet()){
1660 2 : String ephemTab=MeasComet(*(mFrame_p.comet())).getTablePath();
1661 2 : outRecord.define("ephemeristable",ephemTab);
1662 2 : }
1663 48 : return retval;
1664 : }
1665 21 : Bool FTMachine::recoverMovingSourceState(String& error, const RecordInterface& inRecord){
1666 21 : Bool retval=True;
1667 21 : inRecord.get("fixmovingsource", fixMovingSource_p);
1668 21 : { const Record rec=inRecord.asRecord("firstmovingdir_rec");
1669 21 : MeasureHolder mh;
1670 21 : if(!mh.fromRecord(error, rec))
1671 0 : return false;
1672 21 : firstMovingDir_p=mh.asMDirection();
1673 21 : }
1674 21 : { const Record rec=inRecord.asRecord("movingdir_rec");
1675 21 : MeasureHolder mh;
1676 21 : if(!mh.fromRecord(error, rec))
1677 0 : return false;
1678 21 : movingDir_p=mh.asMDirection();
1679 21 : }
1680 21 : { const Record rec=inRecord.asRecord("mlocation_rec");
1681 21 : MeasureHolder mh;
1682 21 : if(!mh.fromRecord(error, rec))
1683 0 : return false;
1684 21 : mLocation_p=mh.asMPosition();
1685 21 : }
1686 21 : SpectralCoordinate *tmpSpec=SpectralCoordinate::restore(inRecord, "spectralcoord");
1687 21 : if(tmpSpec){
1688 21 : spectralCoord_p=*tmpSpec;
1689 21 : delete tmpSpec;
1690 : }
1691 21 : visPolMap_p.resize();
1692 21 : if(inRecord.isDefined("ephemeristable")){
1693 7 : String ephemtab;
1694 7 : inRecord.get("ephemeristable", ephemtab);
1695 7 : MeasComet laComet;
1696 7 : if(Table::isReadable(ephemtab, False)){
1697 7 : Table laTable(ephemtab);
1698 7 : Path leSentier(ephemtab);
1699 7 : laComet=MeasComet(laTable, leSentier.absoluteName());
1700 7 : }
1701 : else{
1702 0 : laComet= MeasComet(ephemtab);
1703 : }
1704 7 : if(!mFrame_p.comet())
1705 7 : mFrame_p.set(laComet);
1706 : else
1707 0 : mFrame_p.resetComet(laComet);
1708 7 : }
1709 :
1710 21 : return retval;
1711 : }
1712 :
1713 :
1714 601329 : void FTMachine::getImagingWeight(Matrix<Float>& imwgt, const vi::VisBuffer2& vb){
1715 : //cerr << "BRIGGSweightor " << briggsWeightor_p.null() << " or " << !briggsWeoght_p << endl;
1716 601329 : if(briggsWeightor_p.null()){
1717 588587 : imwgt=vb.imagingWeight();
1718 : }
1719 : else{
1720 12742 : briggsWeightor_p->weightUniform(imwgt, vb);
1721 : }
1722 :
1723 601329 : }
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 3469 : Bool FTMachine::setFrameValidity(Bool validFrame){
1799 :
1800 3469 : freqFrameValid_p=validFrame;
1801 3469 : return true;
1802 : }
1803 :
1804 :
1805 5534 : Vector<Int> FTMachine::channelMap(const vi::VisBuffer2& vb){
1806 5534 : matchChannel(vb);
1807 5534 : return chanMap;
1808 : }
1809 1040589 : Bool FTMachine::matchChannel(const vi::VisBuffer2& vb){
1810 : //Int spw=vb.spectralWindows()[0];
1811 1040589 : nvischan = vb.nChannels();
1812 1040589 : chanMap.resize(nvischan);
1813 1040589 : chanMap.set(-1);
1814 1040589 : Vector<Double> lsrFreq(0);
1815 :
1816 : //cerr << "doConve " << spw << " " << doConversion_p[spw] << " freqframeval " << freqFrameValid_p << endl;
1817 : //cerr <<"valid frame " << freqFrameValid_p << " polmap "<< polMap << endl;
1818 : //cerr << "spectral coord system " << spectralCoord_p.frequencySystem(False) << endl;
1819 1040589 : if (freqFrameValid_p &&spectralCoord_p.frequencySystem(False)!=MFrequency::REST ) {
1820 821976 : lsrFreq=vb.getFrequencies(0,MFrequency::LSRK);
1821 : }
1822 : else {
1823 218613 : lsrFreq=vb.getFrequencies(0);
1824 : }
1825 1040589 : if (spectralCoord_p.frequencySystem(False)==MFrequency::REST && fixMovingSource_p) {
1826 2393 : if (lastMSId_p != vb.msId()) {
1827 0 : romscol_p=new MSColumns(vb.ms());
1828 : //if ms changed ...reset ephem table
1829 0 : if (upcase(movingDir_p.getRefString()).contains("APP")) {
1830 0 : MeasComet mcomet(Path((romscol_p->field()).ephemPath(vb.fieldId()(0))).absoluteName());
1831 0 : mFrame_p.resetComet(mcomet);
1832 0 : }
1833 : }
1834 4786 : MEpoch e0 (Quantity(vb.time()(0), "s"));
1835 2393 : mFrame_p.epoch() ? mFrame_p.resetEpoch(e0)
1836 0 : : mFrame_p.set(e0);
1837 :
1838 2393 : const auto ephemDirection = vbutil_p->getEphemDir(vb, phaseCenterTime_p);
1839 2393 : mFrame_p.direction() ? mFrame_p.resetDirection(ephemDirection)
1840 0 : : mFrame_p.set(ephemDirection);
1841 :
1842 2393 : shiftFreqToSource(lsrFreq);
1843 2393 : }
1844 : //cerr << "lsrFreq " << lsrFreq.shape() << " nvischan " << nvischan << endl;
1845 : // if(doConversion_p.nelements() < uInt(spw+1))
1846 : // doConversion_p.resize(spw+1, true);
1847 : // doConversion_p[spw]=freqFrameValid_p;
1848 :
1849 1040589 : if(lsrFreq.nelements() ==0){
1850 0 : matchPol(vb);
1851 0 : return false;
1852 : }
1853 1040589 : lsrFreq_p.resize(lsrFreq.nelements());
1854 1040589 : lsrFreq_p=lsrFreq;
1855 1040589 : Vector<Double> c(1);
1856 1040589 : c=0.0;
1857 1040589 : Vector<Double> f(1);
1858 1040589 : Int nFound=0;
1859 :
1860 : Double minFreq;
1861 : Double maxFreq;
1862 1040589 : spectralCoord_p.toWorld(minFreq, Double(0));
1863 1040589 : spectralCoord_p.toWorld(maxFreq, Double(nchan));
1864 1040589 : if(maxFreq < minFreq){
1865 64424 : f(0)=minFreq;
1866 64424 : minFreq=maxFreq;
1867 64424 : maxFreq=f(0);
1868 : }
1869 :
1870 :
1871 : //cout.precision(10);
1872 123133661 : for (Int chan=0;chan<nvischan;chan++) {
1873 122093072 : f(0)=lsrFreq[chan];
1874 122093072 : if(spectralCoord_p.toPixel(c, f)) {
1875 122093072 : Int pixel=Int(floor(c(0)+0.5)); // round to chan freq at chan center
1876 : //cerr << " chan " << chan << " f " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1877 : /////////////
1878 : //c(0)=pixel;
1879 : //spectralCoord_p.toWorld(f, c);
1880 : //cerr << "f1 " << f(0) << " pixel "<< c(0) << " " << pixel << endl;
1881 : ////////////////
1882 122093072 : if(pixel>-1&&pixel<nchan) {
1883 56576142 : chanMap(chan)=pixel;
1884 56576142 : nFound++;
1885 56576142 : if(nvischan>1&&(chan==0||chan==nvischan-1)) {
1886 : /*logIO() << LogIO::DEBUGGING
1887 : << "Selected visibility channel : " << chan+1
1888 : << " has frequency "
1889 : << MFrequency(Quantity(f(0), "Hz")).get("GHz").getValue()
1890 : << " GHz and maps to image pixel " << pixel+1 << LogIO::POST;
1891 : */
1892 : }
1893 : }
1894 : else{
1895 :
1896 65516930 : if(nvischan > 1){
1897 65516930 : Double fwidth=lsrFreq[1]-lsrFreq[0];
1898 65516930 : Double limit=0;
1899 : //Double where=c(0)*fabs(spectralCoord_p.increment()(0));
1900 65516930 : if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::linear)
1901 55770048 : limit=2;
1902 9746882 : else if( freqInterpMethod_p==InterpolateArray1D<Double,Complex>::cubic || freqInterpMethod_p==InterpolateArray1D<Double,Complex>::spline)
1903 0 : limit=4;
1904 : //cerr<< "where " << where << " pixel " << pixel << " fwidth " << fwidth << endl;
1905 : /*
1906 : if(((pixel<0) && (where >= (0-limit*fabs(fwidth)))) )
1907 : chanMap(chan)=-2;
1908 : if((pixel>=nchan) ) {
1909 : where=f(0);
1910 : Double fend;
1911 : spectralCoord_p.toWorld(fend, Double(nchan));
1912 : if( ( (fwidth >0) &&where < (fend+limit*fwidth)) || ( (fwidth <0) &&where > (fend+limit*fwidth)) )
1913 : chanMap(chan)=-2;
1914 : }
1915 : */
1916 :
1917 65516930 : if((f(0) < (maxFreq + limit*fabs(fwidth))) && (f(0) >(maxFreq-0.5*fabs(fwidth)))){
1918 308445 : chanMap(chan)=-2;
1919 : }
1920 65516930 : if((f(0) < minFreq+0.5*fabs(fwidth)) && (f(0) > (minFreq-limit*fabs(fwidth)))){
1921 128373 : chanMap(chan)=-2;
1922 : }
1923 : }
1924 :
1925 :
1926 : }
1927 : }
1928 : }
1929 : //cerr << "chanmap " << chanMap << endl;
1930 : /* if(multiChanMap_p.nelements() < uInt(spw+1))
1931 : multiChanMap_p.resize(spw+1);
1932 : multiChanMap_p[spw].resize();
1933 : multiChanMap_p[spw]=chanMap;
1934 : */
1935 :
1936 1040589 : if(nFound==0) {
1937 : /*
1938 : logIO() << "Visibility channels in spw " << spw+1
1939 : << " of ms " << vb.msId() << " is not being used "
1940 : << LogIO::WARN << LogIO::POST;
1941 : */
1942 19414 : matchPol(vb); ///sometimes the polmap is needed even if chanmap failed
1943 19414 : return false;
1944 : }
1945 :
1946 1021175 : return matchPol(vb);
1947 :
1948 :
1949 :
1950 1040589 : }
1951 :
1952 1040589 : Bool FTMachine::matchPol(const vi::VisBuffer2& vb){
1953 1040589 : Vector<Stokes::StokesTypes> visPolMap(vb.getCorrelationTypesSelected());
1954 1040589 : if((polMap.nelements() > 0) &&(visPolMap.nelements() == visPolMap_p.nelements()) &&allEQ(visPolMap, visPolMap_p))
1955 1036788 : return True;
1956 3801 : Int stokesIndex=image->coordinates().findCoordinate(Coordinate::STOKES);
1957 3801 : AlwaysAssert(stokesIndex>-1, AipsError);
1958 3801 : StokesCoordinate stokesCoord=image->coordinates().stokesCoordinate(stokesIndex);
1959 :
1960 :
1961 3801 : visPolMap_p.resize();
1962 3801 : visPolMap_p=visPolMap;
1963 3801 : nvispol=visPolMap.nelements();
1964 3801 : AlwaysAssert(nvispol>0, AipsError);
1965 3801 : polMap.resize(nvispol);
1966 3801 : polMap=-1;
1967 3801 : Int pol=0;
1968 3801 : Bool found=false;
1969 : // First we try matching Stokes in the visibilities to
1970 : // Stokes in the image that we are gridding into.
1971 12105 : for (pol=0;pol<nvispol;pol++) {
1972 8304 : Int p=0;
1973 8304 : if(stokesCoord.toPixel(p, Stokes::type(visPolMap(pol)))) {
1974 840 : AlwaysAssert(p<npol, AipsError);
1975 840 : polMap(pol)=p;
1976 840 : found=true;
1977 : }
1978 : }
1979 : // If this fails then perhaps we were looking to grid I
1980 : // directly. If so then we need to check that the parallel
1981 : // hands are present in the visibilities.
1982 3801 : if(!found) {
1983 3508 : Int p=0;
1984 3508 : if(stokesCoord.toPixel(p, Stokes::I)) {
1985 3470 : polMap=-1;
1986 3470 : if(vb.polarizationFrame()==MSIter::Linear) {
1987 256 : p=0;
1988 752 : for (pol=0;pol<nvispol;pol++) {
1989 496 : if(Stokes::type(visPolMap(pol))==Stokes::XX)
1990 256 : {polMap(pol)=0;p++;found=true;};
1991 496 : if(Stokes::type(visPolMap(pol))==Stokes::YY)
1992 240 : {polMap(pol)=0;p++;found=true;};
1993 : }
1994 : }
1995 : else {
1996 3214 : p=0;
1997 9898 : for (pol=0;pol<nvispol;pol++) {
1998 6684 : if(Stokes::type(visPolMap(pol))==Stokes::LL)
1999 3214 : {polMap(pol)=0;p++;found=true;};
2000 6684 : if(Stokes::type(visPolMap(pol))==Stokes::RR)
2001 3214 : {polMap(pol)=0;p++;found=true;};
2002 : }
2003 : }
2004 3470 : if(!found) {
2005 : logIO() << "Cannot find polarization map: visibility polarizations = "
2006 0 : << visPolMap << LogIO::EXCEPTION;
2007 : }
2008 : else {
2009 :
2010 : //logIO() << LogIO::DEBUGGING << "Transforming I only" << LogIO::POST;
2011 : }
2012 : };
2013 : }
2014 3801 : return True;
2015 1040589 : }
2016 :
2017 5557 : Vector<String> FTMachine::cleanupTempFiles(const String& mess){
2018 5557 : briggsWeightor_p=nullptr;
2019 5586 : for(uInt k=0; k < tempFileNames_p.nelements(); ++k){
2020 29 : if(Table::isReadable(tempFileNames_p[k])){
2021 29 : if(mess.size()==0){
2022 : try{
2023 1 : Table::deleteTable(tempFileNames_p[k]);
2024 : }
2025 0 : catch(AipsError &x){
2026 0 : logIO() << LogOrigin("FTMachine", "cleanupTempFiles") << LogIO::NORMAL;
2027 0 : logIO() << LogIO::WARN<< "YOU may have to delete the temporary file " << tempFileNames_p[k] << " because " << x.getMesg() << LogIO::POST;
2028 :
2029 0 : }
2030 : }
2031 : else{
2032 28 : logIO() << LogOrigin("FTMachine", "cleanupTempFiles") << LogIO::NORMAL;
2033 28 : logIO() << "YOU have to delete the temporary file " << tempFileNames_p[k] << " because " << mess << LogIO::DEBUG1 << LogIO::POST;
2034 : }
2035 : }
2036 : }
2037 5557 : return tempFileNames_p;
2038 : }
2039 857961 : void FTMachine::gridOk(Int convSupport){
2040 :
2041 857961 : if (nx <= 2*convSupport) {
2042 1 : logIO_p
2043 : << "number of pixels "
2044 : << nx << " on x axis is smaller that the gridding support "
2045 : << 2*convSupport << " Please use a larger value "
2046 1 : << LogIO::EXCEPTION;
2047 : }
2048 :
2049 857960 : if (ny <= 2*convSupport) {
2050 0 : logIO_p
2051 : << "number of pixels "
2052 : << ny << " on y axis is smaller that the gridding support "
2053 : << 2*convSupport << " Please use a larger value "
2054 0 : << LogIO::EXCEPTION;
2055 : }
2056 :
2057 857960 : }
2058 :
2059 0 : void FTMachine::setLocation(const MPosition& loc){
2060 :
2061 0 : mLocation_p=loc;
2062 :
2063 0 : }
2064 :
2065 0 : MPosition& FTMachine::getLocation(){
2066 :
2067 0 : return mLocation_p;
2068 : }
2069 :
2070 :
2071 51 : void FTMachine::setMovingSource(const String& sname, const String& ephtab){
2072 51 : String sourcename=sname;
2073 51 : String ephemtab=ephtab;
2074 : //if a table is given as sourcename...assume ephemerides
2075 51 : if(Table::isReadable(sourcename, False)){
2076 27 : sourcename="COMET";
2077 27 : ephemtab=sname;
2078 27 : ephemTableName_p = sname;
2079 : }
2080 : ///Special case
2081 51 : if(upcase(sourcename)=="TRACKFIELD"){
2082 : //if(name().contains("MosaicFT"))
2083 : // throw(AipsError("Cannot use field phasecenter to track moving source in a Mosaic"));
2084 20 : fixMovingSource_p=True;
2085 20 : movingDir_p=MDirection(Quantity(0.0,"deg"), Quantity(90.0, "deg"));
2086 20 : movingDir_p.setRefString("APP");
2087 : ///Setting it to APP with movingDir_p==True => should use the phasecenter to track
2088 : ///Discussion in CAS-9004 where users are too lazy to give an ephemtable.
2089 20 : return;
2090 : }
2091 :
2092 : MDirection::Types refType;
2093 31 : Bool isValid = MDirection::getType(refType, sourcename);
2094 31 : if(!isValid)
2095 0 : throw(AipsError("Cannot recognize moving source "+sourcename));
2096 31 : if(refType < MDirection::N_Types || refType > MDirection:: N_Planets )
2097 0 : throw(AipsError(sourcename+" is not type of source we can track"));
2098 31 : if(refType==MDirection::COMET){
2099 27 : MeasComet laComet;
2100 27 : if(Table::isReadable(ephemtab, False)){
2101 27 : Table laTable(ephemtab);
2102 27 : Path leSentier(ephemtab);
2103 27 : laComet=MeasComet(laTable, leSentier.absoluteName());
2104 27 : }
2105 : else{
2106 0 : laComet= MeasComet(ephemtab);
2107 : }
2108 27 : if(!mFrame_p.comet())
2109 16 : mFrame_p.set(laComet);
2110 : else
2111 11 : mFrame_p.resetComet(laComet);
2112 27 : }
2113 31 : fixMovingSource_p=true;
2114 31 : movingDir_p=MDirection(Quantity(0.0,"deg"), Quantity(90.0, "deg"));
2115 31 : movingDir_p.setRefString(sourcename);
2116 : // cerr << "movingdir " << movingDir_p.toString() << endl;
2117 71 : }
2118 :
2119 :
2120 0 : void FTMachine::setMovingSource(const MDirection& mdir){
2121 :
2122 0 : fixMovingSource_p=true;
2123 0 : movingDir_p=mdir;
2124 :
2125 0 : }
2126 :
2127 3449 : void FTMachine::setFreqInterpolation(const String& method){
2128 :
2129 3449 : String meth=method;
2130 3449 : meth.downcase();
2131 3449 : if(meth.contains("linear")){
2132 2268 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::linear;
2133 : }
2134 1181 : else if(meth.contains("splin")){
2135 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::spline;
2136 : }
2137 1181 : else if(meth.contains("cub")){
2138 0 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::cubic;
2139 : }
2140 : else{
2141 1181 : freqInterpMethod_p=InterpolateArray1D<Double,Complex>::nearestNeighbour;
2142 : }
2143 :
2144 3449 : }
2145 21 : void FTMachine::setFreqInterpolation(const InterpolateArray1D<Double,Complex>::InterpolationMethod type){
2146 21 : freqInterpMethod_p=type;
2147 21 : }
2148 :
2149 : // helper function to swap the y and z axes of a Cube
2150 339966 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Complex>& in)
2151 : {
2152 339966 : IPosition inShape=in.shape();
2153 339966 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2154 : //resize breaks references...so out better have the right shape
2155 : //if references is not to be broken
2156 339966 : if(out.nelements()==0)
2157 339966 : out.resize(nxx,nyy,nzz);
2158 : Bool deleteIn,deleteOut;
2159 339966 : const Complex* pin = in.getStorage(deleteIn);
2160 339966 : Complex* pout = out.getStorage(deleteOut);
2161 339966 : uInt i=0, zOffset=0;
2162 78634214 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
2163 78294248 : Int yOffset=zOffset;
2164 1570510723 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
2165 4476702225 : for (uInt ix=0; ix<nxx; ++ix){
2166 2984485750 : pout[i++] = pin[ix+yOffset];
2167 : }
2168 : }
2169 : }
2170 339966 : out.putStorage(pout,deleteOut);
2171 339966 : in.freeStorage(pin,deleteIn);
2172 339966 : }
2173 :
2174 53594 : void FTMachine::swapyz(Cube<Complex>& out, const Cube<Bool>& outFlag, const Cube<Complex>& in)
2175 : {
2176 53594 : IPosition inShape=in.shape();
2177 53594 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2178 : //resize breaks references...so out better have the right shape
2179 : //if references is not to be broken
2180 53594 : if(out.nelements()==0)
2181 53594 : out.resize(nxx,nyy,nzz);
2182 : Bool deleteIn,deleteOut, delFlag;
2183 53594 : const Complex* pin = in.getStorage(deleteIn);
2184 53594 : const Bool* poutflag= outFlag.getStorage(delFlag);
2185 53594 : Complex* pout = out.getStorage(deleteOut);
2186 53594 : uInt i=0, zOffset=0;
2187 18850084 : for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
2188 18796490 : Int yOffset=zOffset;
2189 353495373 : for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
2190 1004114249 : for (uInt ix=0; ix<nxx; ++ix){
2191 669415366 : if(!poutflag[i])
2192 538264966 : pout[i] = pin[ix+yOffset];
2193 669415366 : ++i;
2194 : }
2195 : }
2196 : }
2197 53594 : out.putStorage(pout,deleteOut);
2198 53594 : in.freeStorage(pin,deleteIn);
2199 53594 : outFlag.freeStorage(poutflag, delFlag);
2200 53594 : }
2201 :
2202 : // helper function to swap the y and z axes of a Cube
2203 286372 : void FTMachine::swapyz(Cube<Bool>& out, const Cube<Bool>& in)
2204 : {
2205 286372 : IPosition inShape=in.shape();
2206 286372 : uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
2207 286372 : if(out.nelements()==0)
2208 286372 : out.resize(nxx,nyy,nzz);
2209 : Bool deleteIn,deleteOut;
2210 286372 : const Bool* pin = in.getStorage(deleteIn);
2211 286372 : Bool* pout = out.getStorage(deleteOut);
2212 286372 : uInt i=0, zOffset=0;
2213 77788280 : for (uInt iz=0; iz<nzz; iz++, zOffset+=nxx) {
2214 77501908 : Int yOffset=zOffset;
2215 1291907123 : for (uInt iy=0; iy<nyy; iy++, yOffset+=nxx*nzz) {
2216 3643250845 : for (uInt ix=0; ix<nxx; ix++) pout[i++] = pin[ix+yOffset];
2217 : }
2218 : }
2219 286372 : out.putStorage(pout,deleteOut);
2220 286372 : in.freeStorage(pin,deleteIn);
2221 286372 : }
2222 :
2223 137 : void FTMachine::setPointingDirColumn(const String& column){
2224 137 : pointingDirCol_p=column;
2225 137 : pointingDirCol_p.upcase();
2226 137 : if( (pointingDirCol_p != "DIRECTION") &&(pointingDirCol_p != "TARGET") && (pointingDirCol_p != "ENCODER") && (pointingDirCol_p != "POINTING_OFFSET") && (pointingDirCol_p != "SOURCE_OFFSET")){
2227 :
2228 : //basically at this stage you don't know what you're doing...so you get the default
2229 :
2230 0 : pointingDirCol_p="DIRECTION";
2231 :
2232 : }
2233 137 : }
2234 :
2235 0 : String FTMachine::getPointingDirColumnInUse(){
2236 :
2237 0 : return pointingDirCol_p;
2238 :
2239 : }
2240 :
2241 0 : void FTMachine::setSpwChanSelection(const Cube<Int>& spwchansels) {
2242 0 : spwChanSelFlag_p.resize();
2243 0 : spwChanSelFlag_p=spwchansels;
2244 0 : }
2245 :
2246 186 : void FTMachine::setSpwFreqSelection(const Matrix<Double>& spwFreqs)
2247 : {
2248 186 : spwFreqSel_p.assign(spwFreqs);
2249 186 : SynthesisUtils::expandFreqSelection(spwFreqs,expandedSpwFreqSel_p, expandedSpwConjFreqSel_p);
2250 186 : }
2251 :
2252 904752 : void FTMachine::setSpectralFlag(const vi::VisBuffer2& vb, Cube<Bool>& modflagcube){
2253 :
2254 904752 : modflagcube.resize(vb.flagCube().shape());
2255 904752 : modflagcube=vb.flagCube();
2256 904752 : if(!isPseudoI_p){
2257 896813 : ArrayIterator<Bool> from(modflagcube, IPosition(1, 0));
2258 2079565990 : while(!from.pastEnd()){
2259 2078669177 : if(anyTrue(from.array()))
2260 88953118 : from.array().set(True);
2261 2078669177 : from.next();
2262 : }
2263 896813 : }
2264 904752 : uInt nchan = vb.nChannels();
2265 904752 : uInt msid = vb.msId();
2266 904752 : uInt selspw = vb.spectralWindows()[0];
2267 904752 : Bool spwFlagIsSet=( (spwChanSelFlag_p.nelements() > 0) && (spwChanSelFlag_p.shape()(1) > selspw) &&
2268 904752 : (spwChanSelFlag_p.shape()(0) > msid) &&
2269 0 : (spwChanSelFlag_p.shape()(2) >=nchan));
2270 : //cerr << "spwFlagIsSet " << spwFlagIsSet << endl;
2271 101297482 : for (uInt i=0;i<nchan;i++) {
2272 : //Flag those channels that did not get selected...
2273 : //respect the flags from vb if selected or
2274 : //if spwChanSelFlag is wrong shape
2275 100392730 : if ((spwFlagIsSet) && (spwChanSelFlag_p(msid,selspw,i)!=1)) {
2276 0 : modflagcube.xzPlane(i).set(true);
2277 : }
2278 : }
2279 904752 : }
2280 :
2281 789120 : Matrix<Double> FTMachine::negateUV(const vi::VisBuffer2& vb){
2282 789120 : Matrix<Double> uvw(vb.uvw().shape());
2283 252894736 : for (rownr_t i=0;i< vb.nRows() ; ++i) {
2284 756316848 : for (Int idim=0;idim<2; ++idim) uvw(idim,i)=-vb.uvw()(idim, i);
2285 252105616 : uvw(2,i)=vb.uvw()(2,i);
2286 : }
2287 789120 : return uvw;
2288 0 : }
2289 : //-----------------------------------------------------------------------------------------------------------------
2290 : //------------ Vectorized versions of initializeToVis, initializeToSky, finalizeToSky
2291 : //------------ that are called from CubeSkyEquation.
2292 : //------------ They call getImage,getWeightImage, which are implemented in all FTMs
2293 : //------------ Also, Correlation / Stokes conversions and gS/ggS normalizations.
2294 :
2295 :
2296 0 : void FTMachine::setSkyJones(Vector<CountedPtr<casa::refim::SkyJones> >& sj){
2297 0 : sj_p.resize();
2298 0 : sj_p=sj;
2299 0 : cout << "FTM : Set Sky Jones of length " << sj_p.nelements() << " and types ";
2300 0 : for( uInt i=0;i<sj_p.nelements();i++) cout << sj_p[i]->type() << " ";
2301 0 : cout << endl;
2302 0 : }
2303 : // Convert complex correlation planes to float Stokes planes
2304 2763 : void FTMachine::correlationToStokes(ImageInterface<Complex>& compImage,
2305 : ImageInterface<Float>& resImage,
2306 : const Bool dopsf)
2307 : {
2308 : // Convert correlation image to IQUV format
2309 2763 : AlwaysAssert(compImage.shape()[0]==resImage.shape()[0], AipsError);
2310 2763 : AlwaysAssert(compImage.shape()[1]==resImage.shape()[1], AipsError);
2311 2763 : AlwaysAssert(compImage.shape()[3]==resImage.shape()[3], AipsError);
2312 :
2313 2763 : if(dopsf)
2314 : {
2315 : // For the PSF, choose only those stokes planes that have a valid PSF
2316 941 : StokesImageUtil::ToStokesPSF(resImage,compImage);
2317 : }
2318 : else
2319 : {
2320 1822 : StokesImageUtil::To(resImage,compImage);
2321 : }
2322 2763 : };
2323 :
2324 : // Convert float Stokes planes to complex correlation planes
2325 906 : void FTMachine::stokesToCorrelation(ImageInterface<Float>& modelImage,
2326 : ImageInterface<Complex>& compImage)
2327 : {
2328 : /*
2329 : StokesCoordinate stcomp=compImage.coordinates().stokesCoordinate(compImage.coordinates().findCoordinate(Coordinate::STOKES));
2330 : StokesCoordinate stfloat = modelImage.coordinates().stokesCoordinate(modelImage.coordinates().findCoordinate(Coordinate::STOKES));
2331 :
2332 : cout << "Stokes types : complex : " << stcomp.stokes() << " float : " << stfloat.stokes() << endl;
2333 : cout << "Shapes : complex : " << compImage.shape() << " float : " << modelImage.shape() << endl;
2334 : */
2335 :
2336 : //Pol axis need not be same
2337 906 : AlwaysAssert(modelImage.shape()[0]==compImage.shape()[0], AipsError);
2338 906 : AlwaysAssert(modelImage.shape()[1]==compImage.shape()[1], AipsError);
2339 906 : AlwaysAssert(modelImage.shape()[3]==compImage.shape()[3], AipsError);
2340 : // Convert from Stokes to Complex
2341 906 : StokesImageUtil::From(compImage, modelImage);
2342 906 : };
2343 :
2344 : //------------------------------------------------------------------------------------------------------------------
2345 :
2346 0 : void FTMachine::normalizeImage(ImageInterface<Float>& skyImage,
2347 : Matrix<Float>& sumOfWts,
2348 : ImageInterface<Float>& sensitivityImage,
2349 : Bool dopsf, Float pblimit, Int normtype)
2350 : {
2351 :
2352 : //Normalize the sky Image
2353 0 : Int nXX=(skyImage).shape()(0);
2354 0 : Int nYY=(skyImage).shape()(1);
2355 0 : Int npola= (skyImage).shape()(2);
2356 0 : Int nchana= (skyImage).shape()(3);
2357 :
2358 0 : IPosition pcentre(4,nXX/2,nYY/2,0,0);
2359 : // IPosition psource(4,nXX/2+22,nYY/2,0,0);
2360 :
2361 : // storeImg(String("norm_resimage.im") , skyImage);
2362 : // storeImg(String("norm_sensitivity.im"), sensitivityImage);
2363 :
2364 : ///// cout << "FTM::norm : pblimit : " << pblimit << endl;
2365 :
2366 : // Note : This is needed because initial prediction has no info about sumwt.
2367 : // Not a clean solution. // ForSB -- if you see a better way, go for it.
2368 0 : if(sumOfWts.shape() != IPosition(2,npola,nchana))
2369 : {
2370 0 : cout << "Empty SumWt.... resizing and filling with 1.0" << endl;
2371 0 : sumOfWts.resize(IPosition(2,npola,nchana));
2372 0 : sumOfWts=1.0;
2373 : }
2374 :
2375 : // if(dopsf) cout << "*** FTM::normalizeImage : Image Center : " << skyImage.getAt(pcentre) << " and weightImage : " << sensitivityImage.getAt(pcentre) << " SumWt : " << sumOfWts[0,0];
2376 : // else cout << "*** FTM::normalizeImage : Source Loc : " << skyImage.getAt(psource) << " and weightImage : " << sensitivityImage.getAt(psource) << " SumWt : " << sumOfWts[0,0];
2377 :
2378 :
2379 :
2380 0 : IPosition blc(4,nXX, nYY, npola, nchana);
2381 0 : IPosition trc(4, nXX, nYY, npola, nchana);
2382 0 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
2383 : //max weights per plane
2384 0 : for (Int pol=0; pol < npola; ++pol){
2385 0 : for (Int chan=0; chan < nchana ; ++chan){
2386 :
2387 0 : blc(2)=pol; trc(2)=pol;
2388 0 : blc(3)=chan; trc(3)=chan;
2389 0 : Slicer sl(blc, trc, Slicer::endIsLast);
2390 0 : SubImage<Float> subSkyImage(skyImage, sl, false);
2391 0 : SubImage<Float> subSensitivityImage(sensitivityImage, sl, false);
2392 0 : SubImage<Float> subOutput(skyImage, sl, true);
2393 0 : Float sumWt = sumOfWts(pol,chan);
2394 0 : if(sumWt > 0.0){
2395 0 : switch(normtype)
2396 : {
2397 0 : case 0: // only sum Of Weights - FTM only (ForSB)
2398 0 : subOutput.copyData( (LatticeExpr<Float>)
2399 0 : ((dopsf?1.0:-1.0)*subSkyImage/(sumWt)));
2400 0 : break;
2401 :
2402 0 : case 1: // only sensitivityImage Ic/avgPB (ForSB)
2403 0 : subOutput.copyData( (LatticeExpr<Float>)
2404 0 : (iif(subSensitivityImage > (pblimit),
2405 0 : (subSkyImage/(subSensitivityImage)),
2406 : (subSkyImage))));
2407 : // 0.0)));
2408 0 : break;
2409 :
2410 0 : case 2: // sum of Weights and sensitivityImage IGridded/(SoW*avgPB) and PSF --> Id (ForSB)
2411 0 : subOutput.copyData( (LatticeExpr<Float>)
2412 0 : (iif(subSensitivityImage > (pblimit),
2413 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage*sumWt)),
2414 : //((dopsf?1.0:-1.0)*subSkyImage))));
2415 : 0.0)));
2416 0 : break;
2417 :
2418 0 : case 3: // MULTIPLY by the sensitivityImage avgPB
2419 0 : subOutput.copyData( (LatticeExpr<Float>) (subSkyImage * subSensitivityImage) );
2420 0 : break;
2421 :
2422 0 : case 4: // DIVIDE by sqrt of sensitivityImage
2423 0 : subOutput.copyData( (LatticeExpr<Float>)
2424 0 : (iif((subSensitivityImage) > (pblimit),
2425 0 : (subSkyImage/(sqrt(subSensitivityImage))),
2426 : (subSkyImage))));
2427 : //0.0)));
2428 0 : break;
2429 :
2430 0 : case 5: // MULTIPLY by sqrt of sensitivityImage
2431 0 : subOutput.copyData( (LatticeExpr<Float>)
2432 0 : (iif((subSensitivityImage) > (pblimit),
2433 0 : (subSkyImage * (sqrt(subSensitivityImage))),
2434 : (subSkyImage))));
2435 :
2436 0 : break;
2437 :
2438 0 : case 6: // divide by non normalized sensitivity image
2439 : {
2440 0 : Float elpblimit=max( subSensitivityImage).getFloat() * pblimit;
2441 0 : subOutput.copyData( (LatticeExpr<Float>)
2442 0 : (iif(subSensitivityImage > (elpblimit),
2443 0 : ((dopsf?1.0:-1.0)*subSkyImage/(subSensitivityImage)),
2444 : 0.0)));
2445 : }
2446 0 : break;
2447 0 : default:
2448 0 : throw(AipsError("Unrecognized norm-type in FTM::normalizeImage : "+String::toString(normtype)));
2449 : }
2450 : }
2451 : else{
2452 0 : subOutput.set(0.0);
2453 : }
2454 0 : }
2455 : }
2456 :
2457 : //if(dopsf) cout << " Normalized (" << normtype << ") Image Center : " << skyImage.getAt(pcentre) << endl;
2458 : // else cout << " Normalized (" << normtype << ") Source Loc : " << skyImage.getAt(psource) << endl;
2459 :
2460 0 : }
2461 :
2462 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2463 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2464 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2465 : ///// For use with the new framework
2466 : ///// (Sorry about these copies, but need to keep old system working)
2467 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2468 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2469 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
2470 :
2471 : // Vectorized InitializeToVis
2472 750 : void FTMachine::initializeToVisNew(const VisBuffer2& vb,
2473 : CountedPtr<SIImageStore> imstore)
2474 : {
2475 750 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2476 :
2477 750 : Matrix<Float> tempWts;
2478 :
2479 750 : if(!(imstore->forwardGrid()).get())
2480 0 : throw(AipsError("FTMAchine::InitializeToVisNew error imagestore has no valid grid initialized"));
2481 : // Convert from Stokes planes to Correlation planes
2482 750 : LatticeLocker lock1 (*(imstore->model()), FileLocker::Read);
2483 750 : stokesToCorrelation(*(imstore->model()), *(imstore->forwardGrid()));
2484 :
2485 750 : if(vb.polarizationFrame()==MSIter::Linear) {
2486 50 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2487 : StokesImageUtil::LINEAR);
2488 : }
2489 : else {
2490 700 : StokesImageUtil::changeCStokesRep(*(imstore->forwardGrid()),
2491 : StokesImageUtil::CIRCULAR);
2492 : }
2493 :
2494 : //------------------------------------------------------------------------------------
2495 : // Image Mosaic only : Multiply the input model with the Primary Beam
2496 750 : if(sj_p.nelements() >0 ){
2497 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2498 0 : (sj_p(k))->apply(*(imstore->forwardGrid()), *(imstore->forwardGrid()), vb, 0, true);
2499 : }
2500 : }
2501 : //------------------------------------------------------------------------------------
2502 :
2503 : // Call initializeToVis
2504 750 : initializeToVis(*(imstore->forwardGrid()), vb); // Pure virtual
2505 :
2506 750 : };
2507 :
2508 : // Vectorized finalizeToVis is not implemented because it does nothing and is never called.
2509 :
2510 : // Vectorized InitializeToSky
2511 2097 : void FTMachine::initializeToSkyNew(const Bool dopsf,
2512 : const VisBuffer2& vb,
2513 : CountedPtr<SIImageStore> imstore)
2514 :
2515 : {
2516 2097 : AlwaysAssert(imstore->getNTaylorTerms(false)==1, AipsError);
2517 :
2518 : // Make the relevant float grid.
2519 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
2520 2097 : if( dopsf ) { imstore->psf(); } else { imstore->residual(); }
2521 :
2522 : // Initialize the complex grid (i.e. tell FTMachine what array to use internally)
2523 2096 : Matrix<Float> sumWeight;
2524 2096 : if(!(imstore->backwardGrid()).get())
2525 0 : throw(AipsError("FTMAchine::InitializeToSkyNew error imagestore has no valid grid initialized"));
2526 2096 : initializeToSky(*(imstore->backwardGrid()) , sumWeight , vb);
2527 :
2528 2096 : };
2529 :
2530 : // Vectorized finalizeToSky
2531 2088 : void FTMachine::finalizeToSkyNew(Bool dopsf,
2532 : const VisBuffer2& vb,
2533 : CountedPtr<SIImageStore> imstore )
2534 : {
2535 : // Check vector lengths.
2536 2088 : AlwaysAssert( imstore->getNTaylorTerms(false)==1, AipsError);
2537 :
2538 2088 : Matrix<Float> sumWeights;
2539 2088 : finalizeToSky();
2540 :
2541 : //------------------------------------------------------------------------------------
2542 : // Straightforward case. No extra primary beams. No image mosaic
2543 2088 : if(sj_p.nelements() == 0 )
2544 : {
2545 : // cerr << "TYPEID " << typeid( *(imstore->psf())).name() << " " << typeid(typeid( *(imstore->psf())).name()).name() << endl;
2546 2088 : shared_ptr<ImageInterface<Float> > theim=dopsf ? imstore->psf() : imstore->residual();
2547 : //static_cast<decltype(imstore->residual())>(theim)->lock();
2548 2088 : { LatticeLocker lock1 (*theim, FileLocker::Write);
2549 2088 : correlationToStokes( getImage(sumWeights, false) , *theim, dopsf);
2550 2088 : }
2551 2088 : theim->unlock();
2552 :
2553 2088 : if( (useWeightImage() && dopsf) || isSD() ) {
2554 :
2555 207 : LatticeLocker lock1 (*(imstore->weight()), FileLocker::Write);
2556 207 : getWeightImage( *(imstore->weight()) , sumWeights);
2557 207 : imstore->weight()->unlock();
2558 :
2559 : // Fill weight image only once, during PSF generation. Remember.... it is normalized only once
2560 : // during PSF generation.
2561 207 : }
2562 :
2563 : // Take sumWeights from corrToStokes here....
2564 2088 : LatticeLocker lock1 (*(imstore->sumwt()), FileLocker::Write);
2565 2088 : Bool donesumwt=(max(imstore->sumwt()->get()) > 0.0);
2566 2088 : if(!donesumwt){
2567 1478 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
2568 739 : CoordinateSystem incoord=image->coordinates();
2569 739 : CoordinateSystem outcoord=imstore->sumwt()->coordinates();
2570 739 : StokesImageUtil::ToStokesSumWt(sumWeightStokes, sumWeights, outcoord, incoord);
2571 :
2572 :
2573 739 : Array<Float> sumWtArr(IPosition(4,1,1,sumWeights.shape()[0], sumWeights.shape()[1]));
2574 :
2575 739 : IPosition blc(4, 0, 0, 0, 0);
2576 739 : IPosition trc(4, 0, 0, sumWeightStokes.shape()[0]-1, sumWeightStokes.shape()[1]-1);
2577 739 : sumWtArr(blc, trc).reform(sumWeightStokes.shape())=sumWeightStokes;
2578 :
2579 : //StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
2580 739 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeightStokes.shape()[0] ) &&
2581 : ((imstore->sumwt())->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
2582 :
2583 739 : (imstore->sumwt())->put( sumWeightStokes.reform((imstore->sumwt())->shape()) );
2584 739 : }
2585 2088 : imstore->sumwt()->unlock();
2586 :
2587 2088 : }
2588 : //------------------------------------------------------------------------------------
2589 : // Image Mosaic only : Multiply the residual, and weight image by the PB.
2590 : else
2591 : {
2592 :
2593 : // Take the FT of the gridded values. Writes into backwardGrid().
2594 0 : getImage(sumWeights, false);
2595 :
2596 : // Multiply complex image grid by PB.
2597 0 : if( !dopsf )
2598 : {
2599 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2600 0 : (sj_p(k))->apply(*(imstore->backwardGrid()), *(imstore->backwardGrid()), vb, 0, true);
2601 : }
2602 : }
2603 :
2604 : // Convert from correlation to Stokes onto a new temporary grid.
2605 0 : SubImage<Float> targetImage( ( dopsf ? *(imstore->psf()) : *(imstore->residual()) ) , true);
2606 0 : TempImage<Float> temp( targetImage.shape(), targetImage.coordinates() );
2607 0 : correlationToStokes( *(imstore->backwardGrid()), temp, false);
2608 :
2609 : // Add the temporary Stokes image to the residual or PSF, whichever is being made.
2610 0 : LatticeExpr<Float> addToRes( targetImage + temp );
2611 0 : targetImage.copyData(addToRes);
2612 :
2613 : // Now, do the same with the weight image and sumwt ( only on the first pass )
2614 0 : if( dopsf )
2615 : {
2616 0 : SubImage<Float> weightImage( *(imstore->weight()) , true);
2617 0 : TempImage<Float> temp(weightImage.shape(), weightImage.coordinates());
2618 0 : getWeightImage(temp, sumWeights);
2619 :
2620 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2621 0 : (sj_p(k))->applySquare(temp,temp, vb, -1);
2622 : }
2623 :
2624 0 : LatticeExpr<Float> addToWgt( weightImage + temp );
2625 0 : weightImage.copyData(addToWgt);
2626 :
2627 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeights.shape()[0] ) &&
2628 : ((imstore->sumwt())->shape()[3] == sumWeights.shape()[1] ) , AipsError );
2629 :
2630 0 : SubImage<Float> sumwtImage( *(imstore->sumwt()) , true);
2631 0 : TempImage<Float> temp2(sumwtImage.shape(), sumwtImage.coordinates());
2632 0 : temp2.put( sumWeights.reform(sumwtImage.shape()) );
2633 0 : LatticeExpr<Float> addToWgt2( sumwtImage + temp2 );
2634 0 : sumwtImage.copyData(addToWgt2);
2635 :
2636 : //cout << "In finalizeGridCoreMos : sumwt : " << sumwtImage.get() << endl;
2637 :
2638 0 : }
2639 :
2640 0 : }
2641 : //------------------------------------------------------------------------------------
2642 :
2643 :
2644 :
2645 4176 : return;
2646 2088 : };
2647 :
2648 :
2649 : /////------------------------------------------------
2650 0 : void FTMachine::finalizeToWeightImage(const VisBuffer2& vb,
2651 : CountedPtr<SIImageStore> imstore )
2652 : {
2653 : // Check vector lengths.
2654 0 : AlwaysAssert( imstore->getNTaylorTerms(false)==1, AipsError);
2655 :
2656 0 : Matrix<Float> sumWeights;
2657 :
2658 : //------------------------------------------------------------------------------------
2659 : // Straightforward case. No extra primary beams. No image mosaic
2660 0 : if(sj_p.nelements() == 0 )
2661 : {
2662 :
2663 :
2664 0 : if( useWeightImage() ) {
2665 : //if( name().contains("Mosaic") ){
2666 : {
2667 0 : finalizeToSky();
2668 : }
2669 0 : LatticeLocker lock1 (*(imstore->weight()), FileLocker::Write);
2670 0 : getWeightImage( *(imstore->weight()) , sumWeights);
2671 0 : imstore->weight()->unlock();
2672 :
2673 : // Fill weight image only once, during PSF generation. Remember.... it is normalized only once
2674 : // during PSF generation.
2675 0 : }
2676 0 : if(sumWeights.nelements() >0){
2677 : // Take sumWeights from corrToStokes here....
2678 0 : LatticeLocker lock1 (*(imstore->sumwt()), FileLocker::Write);
2679 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
2680 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
2681 :
2682 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeightStokes.shape()[0] ) &&
2683 : ((imstore->sumwt())->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
2684 :
2685 0 : (imstore->sumwt())->put( sumWeightStokes.reform((imstore->sumwt())->shape()) );
2686 0 : imstore->sumwt()->unlock();
2687 0 : }
2688 :
2689 : }
2690 : //------------------------------------------------------------------------------------
2691 : // Image Mosaic only : Multiply the residual, and weight image by the PB.
2692 : else
2693 : {
2694 :
2695 : // Now, do the same with the weight image and sumwt ( only on the first pass )
2696 : {
2697 0 : SubImage<Float> weightImage( *(imstore->weight()) , true);
2698 0 : TempImage<Float> temp(weightImage.shape(), weightImage.coordinates());
2699 0 : getWeightImage(temp, sumWeights);
2700 :
2701 0 : for (uInt k=0; k < sj_p.nelements(); ++k){
2702 0 : (sj_p(k))->applySquare(temp,temp, vb, -1);
2703 : }
2704 :
2705 0 : LatticeExpr<Float> addToWgt( weightImage + temp );
2706 0 : weightImage.copyData(addToWgt);
2707 :
2708 0 : AlwaysAssert( ( (imstore->sumwt())->shape()[2] == sumWeights.shape()[0] ) &&
2709 : ((imstore->sumwt())->shape()[3] == sumWeights.shape()[1] ) , AipsError );
2710 :
2711 0 : SubImage<Float> sumwtImage( *(imstore->sumwt()) , true);
2712 0 : TempImage<Float> temp2(sumwtImage.shape(), sumwtImage.coordinates());
2713 0 : temp2.put( sumWeights.reform(sumwtImage.shape()) );
2714 0 : LatticeExpr<Float> addToWgt2( sumwtImage + temp2 );
2715 0 : sumwtImage.copyData(addToWgt2);
2716 :
2717 : //cout << "In finalizeGridCoreMos : sumwt : " << sumwtImage.get() << endl;
2718 :
2719 0 : }
2720 :
2721 : }
2722 : //------------------------------------------------------------------------------------
2723 :
2724 :
2725 :
2726 0 : return;
2727 0 : };
2728 :
2729 :
2730 :
2731 :
2732 : /////-----------------------------------------------
2733 0 : Bool FTMachine::changedSkyJonesLogic(const vi::VisBuffer2& vb, Bool& firstRow, Bool& internalRow)
2734 : {
2735 0 : firstRow=false;
2736 0 : internalRow=false;
2737 :
2738 0 : if( sj_p.nelements()==0 )
2739 0 : {throw(AipsError("Internal Error : Checking changedSkyJones, but it is not yet set."));}
2740 :
2741 0 : CountedPtr<SkyJones> ej = sj_p[0];
2742 0 : if(ej.null())
2743 0 : return false;
2744 0 : if(ej->changed(vb,0))
2745 0 : firstRow=true;
2746 0 : Int row2temp=0;
2747 0 : if(ej->changedBuffer(vb,0,row2temp)) {
2748 0 : internalRow=true;
2749 : }
2750 0 : return (firstRow || internalRow) ;
2751 0 : }
2752 :
2753 0 : std::shared_ptr<std::complex<double>> FTMachine::getGridPtr(size_t& size) const
2754 : {
2755 0 : size = 0;
2756 0 : return std::shared_ptr<std::complex<double>>();
2757 : }
2758 :
2759 0 : std::shared_ptr<double> FTMachine::getSumWeightsPtr(size_t& size) const
2760 : {
2761 0 : size = 0;
2762 0 : return std::shared_ptr<double>();
2763 : }
2764 :
2765 0 : void FTMachine::setCFCache(CountedPtr<CFCache>& /*cfc*/, const Bool /*loadCFC*/)
2766 : {
2767 0 : throw(AipsError("FTMachine::setCFCache() directly called!"));
2768 : }
2769 :
2770 :
2771 :
2772 49464 : 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){
2773 : /* Vector<Int> ord(36);
2774 : ord(0)=14;
2775 : ord(1)=15;
2776 : ord(2)=20;
2777 : ord(3)=21;ord(4)=13;
2778 : ord(5)=16;ord(6)=19;ord(7)=22;ord(8)=8;ord(9)=9;
2779 : ord(10)=26;ord(11)=27;ord(12)=25;ord(13)=28;ord(14)=7;
2780 : ord(15)=10;ord(16)=32;ord(17)=33;ord(18)=2;ord(19)=3;
2781 : ord(20)=18;ord(21)=23;ord(22)=12;ord(23)=17;ord(24)=1;
2782 : ord(25)=4;ord(26)=6;ord(27)=11;ord(28)=24;ord(29)=29;
2783 : ord(30)=31;ord(31)=34;ord(32)=0;ord(33)=5;ord(34)=30;
2784 : ord(35)=35;
2785 : */
2786 : /*
2787 : Int ix= (icounter+1)-((icounter)/ixsub)*ixsub;
2788 : Int iy=(icounter)/ixsub+1;
2789 : y0=(nyp/iysub)*(iy-1)+1+miny;
2790 : nysub=nyp/iysub;
2791 : if( iy == iysub) {
2792 : nysub=nyp-(nyp/iysub)*(iy-1);
2793 : }
2794 : x0=(nxp/ixsub)*(ix-1)+1+minx;
2795 : nxsub=nxp/ixsub;
2796 : if( ix == ixsub){
2797 : nxsub=nxp-(nxp/ixsub)*(ix-1);
2798 : }
2799 : */
2800 :
2801 :
2802 : {
2803 49464 : Int elrow=icounter/ixsub;
2804 49464 : Int elcol=(icounter-elrow*ixsub);
2805 : //cerr << "row "<< elrow << " col " << elcol << endl;
2806 : //nxsub=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0 + 0.5));
2807 49464 : Float factor=0;
2808 49464 : if(ixsub > 1){
2809 238896 : for (Int k=0; k < ixsub/2; ++k)
2810 189432 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
2811 : //factor= linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
2812 49464 : factor *= 2.0;
2813 49464 : if(linear)
2814 49464 : nxsub=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2815 : else
2816 : //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));
2817 0 : nxsub=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
2818 : }
2819 : else{
2820 0 : nxsub=nxp;
2821 : }
2822 : //cerr << nxp << " col " << elcol << " nxsub " << nxsub << endl;
2823 49464 : x0=minx;
2824 49464 : elcol-=1;
2825 214164 : while(elcol >= 0){
2826 : //x0+=Int(floor(((ceil(fabs(float(2*elcol+1-ixsub)/2.0))-1.0)*5 +1)*nxp/36.0+0.5));
2827 :
2828 164700 : if(linear)
2829 164700 : x0+=Int(floor((ceil(fabs(float(2*elcol+1-ixsub)/2.0))/factor)*nxp + 0.5));
2830 : else
2831 : //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));
2832 0 : x0+=Int(floor((sqrt(ceil(fabs(float(2*elcol+1-ixsub)/2.0)))/factor)*nxp + 0.5));
2833 164700 : elcol-=1;
2834 : }
2835 49464 : factor=0;
2836 49464 : if(iysub >1){
2837 238896 : for (Int k=0; k < iysub/2; ++k)
2838 : //factor=linear ? factor+(k+1): factor+(k+1)*(k+1)*(k+1);
2839 189432 : factor= linear ? factor+(k+1): factor+sqrt(Float(k+1));
2840 49464 : factor *= 2.0;
2841 : //nysub=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
2842 49464 : if(linear)
2843 49464 : nysub=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2844 : else
2845 0 : nysub=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
2846 : }
2847 : else{
2848 0 : nysub=nyp;
2849 : }
2850 : //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));
2851 49464 : y0=miny;
2852 49464 : elrow-=1;
2853 :
2854 214164 : while(elrow >=0){
2855 : //y0+=Int(floor(((ceil(fabs(float(2*elrow+1-iysub)/2.0))-1.0)*5 +1)*nyp/36.0+0.5));
2856 164700 : if(linear)
2857 164700 : y0+=Int(floor((ceil(fabs(float(2*elrow+1-iysub)/2.0))/factor)*nyp + 0.5));
2858 : else
2859 0 : y0+=Int(floor((sqrt(ceil(fabs(float(2*elrow+1-iysub)/2.0)))/factor)*nyp + 0.5));
2860 : //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));
2861 164700 : elrow-=1;
2862 : }
2863 : }
2864 49464 : if( (y0+nysub) >= nyp){
2865 7236 : nysub = nyp-y0-1;
2866 : }
2867 49464 : if( (x0+nxsub) >= nxp){
2868 7236 : nxsub = nxp-x0-1;
2869 : }
2870 49464 : y0+=1;
2871 49464 : x0+=1;
2872 : //cerr << icounter << " x0, y0 " << x0 << " " << y0 << " ixsub, iysub " << nxsub << " " << nysub << endl;
2873 49464 : if(doneThreadPartition_p < 0)
2874 1431 : doneThreadPartition_p=1;
2875 :
2876 49464 : }
2877 :
2878 111036 : void FTMachine::tweakGridSector(const Int& nx, const Int& ny, const Int& ixsub, const Int& iysub){
2879 : //if(doneThreadPartition_p)
2880 : // return;
2881 111036 : Vector<Int> x0, y0, nxsub, nysub;
2882 111036 : Vector<Float> xcut(nx/2);
2883 111036 : Vector<Float> ycut(ny/2);
2884 111036 : if(griddedData2.nelements() >0 ){
2885 : //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;
2886 110556 : 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()));
2887 110556 : 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()));
2888 : }
2889 : else{
2890 480 : xcut=amplitude(griddedData(IPosition(4, 0, ny/2-1, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
2891 480 : ycut=amplitude(griddedData(IPosition(4, nx/2-1, 0, 0, 0), IPosition(4, nx/2-1, ny/2-1, 0,0)));
2892 : }
2893 : //cerr << griddedData2.shape() << " " << griddedData.shape() << endl;
2894 111036 : Vector<Float> cumSumX(nx/2, 0);
2895 : //Vector<Float> cumSumX2(nx/2,0);
2896 111036 : cumSumX(0)=xcut(0);
2897 : //cumSumX2(0)=cumSumX(0)*cumSumX(0);
2898 111036 : Float sumX=sum(xcut);
2899 111036 : if(sumX==0.0)
2900 1732 : return;
2901 : //cerr << "sumX " << sumX << endl;
2902 : //sumX *= sumX;
2903 109304 : x0.resize(ixsub);
2904 109304 : x0=nx/2-1;
2905 109304 : nxsub.resize(ixsub);
2906 109304 : nxsub=0;
2907 109304 : x0(0)=0;
2908 109304 : Int counter=1;
2909 14759344 : for (Int k=1; k < nx/2; ++k){
2910 14650040 : cumSumX(k)=cumSumX(k-1)+xcut(k);
2911 : //cumSumX2(k)=cumSumX(k)*cumSumX(k);
2912 14650040 : Float nextEdge=sumX/(Float(ixsub/2)*Float(ixsub/2))*Float(counter)*Float(counter);
2913 14650040 : if(cumSumX(k-1) < nextEdge && cumSumX(k) >= nextEdge){
2914 378231 : x0(counter)=k;
2915 : //cerr << counter << " " << k << " diff " << x0(counter)-x0(counter-1) << endl;
2916 378231 : nxsub(counter-1)=x0(counter)-x0(counter-1);
2917 378231 : ++counter;
2918 : }
2919 : }
2920 :
2921 109304 : x0(ixsub/2)=nx/2-1;
2922 109304 : nxsub(ixsub/2)=nxsub(ixsub/2-1)=nx/2-1-x0(ixsub/2-1);
2923 437216 : for(Int k=ixsub/2+1; k < ixsub; ++k){
2924 327912 : x0(k)=x0(k-1)+ nxsub(ixsub-k);
2925 327912 : nxsub(k)=nxsub(ixsub-k-1);
2926 : }
2927 109304 : nxsub(ixsub-1)+=1;
2928 :
2929 109304 : Vector<Float> cumSumY(ny/2, 0);
2930 : //Vector<Float> cumSumY2(ny/2,0);
2931 109304 : cumSumY(0)=ycut(0);
2932 : //cumSumY2(0)=cumSumY(0)*cumSumY(0);
2933 109304 : Float sumY=sum(ycut);
2934 109304 : if(sumY==0.0)
2935 33 : return;
2936 : //sumY *=sumY;
2937 109271 : y0.resize(iysub);
2938 109271 : y0=ny/2-1;
2939 109271 : nysub.resize(iysub);
2940 109271 : nysub=0;
2941 109271 : y0(0)=0;
2942 109271 : counter=1;
2943 14757364 : for (Int k=1; k < ny/2; ++k){
2944 14648093 : cumSumY(k)=cumSumY(k-1)+ycut(k);
2945 : //cumSumY2(k)=cumSumY(k)*cumSumY(k);
2946 14648093 : Float nextEdge=sumY/(Float(iysub/2)*Float(iysub/2))*Float(counter)*Float(counter);
2947 14648093 : if(cumSumY(k-1) < nextEdge && cumSumY(k) >= nextEdge){
2948 399230 : y0(counter)=k;
2949 399230 : nysub(counter-1)=y0(counter)-y0(counter-1);
2950 399230 : ++counter;
2951 : }
2952 : }
2953 :
2954 109271 : y0(ixsub/2)=ny/2-1;
2955 109271 : nysub(iysub/2)=nysub(iysub/2-1)=ny/2-1-y0(iysub/2-1);
2956 437084 : for(Int k=iysub/2+1; k < iysub; ++k){
2957 327813 : y0(k)=y0(k-1)+ nysub(iysub-k);
2958 327813 : nysub(k)=nysub(iysub-k-1);
2959 : }
2960 109271 : nysub(iysub-1)+=1;
2961 :
2962 109271 : if(anyEQ(nxsub, 0) || anyEQ(nysub, 0))
2963 69895 : return;
2964 : //cerr << " x0 " << x0 << " nxsub " << nxsub << endl;
2965 : //cerr << " y0 " << y0 << " nysub " << nysub << endl;
2966 39376 : x0+=1;
2967 39376 : y0+=1;
2968 39376 : xsect_p.resize(ixsub*iysub);
2969 39376 : ysect_p.resize(ixsub*iysub);
2970 39376 : nxsect_p.resize(ixsub*iysub);
2971 39376 : nysect_p.resize(ixsub*iysub);
2972 354384 : for (Int iy=0; iy < iysub; ++iy){
2973 2835072 : for (Int ix=0; ix< ixsub; ++ix){
2974 :
2975 2520064 : xsect_p(iy*ixsub+ix)=x0[ix];
2976 2520064 : ysect_p(iy*ixsub+ix)=y0[iy];
2977 2520064 : nxsect_p(iy*ixsub+ix)=nxsub[ix];
2978 2520064 : nysect_p(iy*ixsub+ix)=nysub[iy];
2979 : }
2980 : }
2981 :
2982 39376 : ++doneThreadPartition_p;
2983 :
2984 610924 : }
2985 :
2986 :
2987 : /*
2988 : /// Move to individual FTMs............ make it pure virtual.
2989 : Bool FTMachine::useWeightImage()
2990 : {
2991 : if( name() == "GridFT" || name() == "WProjectFT" )
2992 : { return false; }
2993 : else
2994 : { return true; }
2995 : }
2996 : */
2997 :
2998 : }//# namespace refim ends
2999 : }//namespace CASA ends
3000 :
|
