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