Commits

Takeshi Nakazato authored and Ville Suoranta committed 66bed521c76 Merge
Pull request #777: CAS-14266

Merge in CASA/casa6 from CAS-14266 to master * commit '35aa5b675a4d465b976538cdcd112c1f70eb3772': CAS-14266 do not initialize UVWMachine for single-dish imaging

casatools/src/code/synthesis/TransformMachines2/FTMachine.cc

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

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut