Line data Source code
1 : // -*- C++ -*-
2 : //# Imager.cc: Implementation of Imager.h
3 : //# All helper functions of imager moved here for readability
4 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
5 : //# Associated Universities, Inc. Washington DC, USA.
6 : //#
7 : //# This library is free software; you can redistribute it and/or modify it
8 : //# under the terms of the GNU Library General Public License as published by
9 : //# the Free Software Foundation; either version 2 of the License, or (at your
10 : //# option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful, but WITHOUT
13 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15 : //# License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Library General Public License
18 : //# along with this library; if not, write to the Free Software Foundation,
19 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20 : //#
21 : //# Correspondence concerning AIPS++ should be addressed as follows:
22 : //# Internet email: casa-feedback@nrao.edu.
23 : //# Postal address: AIPS++ Project Office
24 : //# National Radio Astronomy Observatory
25 : //# 520 Edgemont Road
26 : //# Charlottesville, VA 22903-2475 USA
27 : //#
28 : //# $Id$
29 :
30 :
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <iostream>
33 : #include <synthesis/MeasurementEquations/Imager.h>
34 : #include <synthesis/MeasurementComponents/EPJones.h>
35 : #include <synthesis/TransformMachines/VisModelData.h>
36 :
37 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
38 :
39 : #include <casacore/casa/Arrays/Matrix.h>
40 : #include <casacore/casa/Arrays/ArrayMath.h>
41 : #include <casacore/casa/Arrays/ArrayLogical.h>
42 :
43 : #include <casacore/casa/Logging.h>
44 : #include <casacore/casa/Logging/LogIO.h>
45 : #include <casacore/casa/Logging/LogMessage.h>
46 :
47 : #include <casacore/casa/OS/DirectoryIterator.h>
48 : #include <casacore/casa/OS/File.h>
49 : #include <casacore/casa/OS/Path.h>
50 :
51 : #include <casacore/casa/OS/HostInfo.h>
52 : #include <casacore/tables/Tables/RefRows.h>
53 : #include <casacore/tables/Tables/Table.h>
54 : #include <casacore/tables/Tables/SetupNewTab.h>
55 : #include <casacore/tables/TaQL/TableParse.h>
56 : #include <casacore/tables/Tables/TableRecord.h>
57 : #include <casacore/tables/Tables/TableDesc.h>
58 : #include <casacore/tables/Tables/TableLock.h>
59 : #include <casacore/tables/TaQL/ExprNode.h>
60 :
61 : #include <casacore/casa/BasicSL/String.h>
62 : #include <casacore/casa/Utilities/Assert.h>
63 : #include <casacore/casa/Utilities/Fallible.h>
64 : #include <casacore/casa/Utilities/CompositeNumber.h>
65 :
66 : #include <casacore/casa/BasicSL/Constants.h>
67 :
68 : #include <casacore/casa/Logging/LogSink.h>
69 : #include <casacore/casa/Logging/LogMessage.h>
70 :
71 : #include <casacore/casa/Arrays/ArrayMath.h>
72 : #include <casacore/casa/Arrays/Slice.h>
73 : #include <casacore/images/Images/ImageExpr.h>
74 : #include <imageanalysis/ImageAnalysis/ImagePolarimetry.h>
75 : #include <imageanalysis/Images/ComponentListImage.h>
76 : #include <casacore/images/Images/ImageBeamSet.h>
77 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
78 : #include <casacore/lattices/LatticeMath/LatticeCleanProgress.h>
79 : #include <msvis/MSVis/VisSet.h>
80 : #include <msvis/MSVis/VisSetUtil.h>
81 : #include <msvis/MSVis/VisImagingWeight.h>
82 : #include <msvis/MSVis/SubMS.h>
83 : // Disabling Imager::correct() (gmoellen 06Nov20)
84 : //#include <synthesis/MeasurementComponents/TimeVarVisJones.h>
85 :
86 : #include <casacore/measures/Measures/Stokes.h>
87 : #include <casacore/casa/Quanta/UnitMap.h>
88 : #include <casacore/casa/Quanta/UnitVal.h>
89 : #include <casacore/casa/Quanta/MVAngle.h>
90 : #include <casacore/measures/Measures/MDirection.h>
91 : #include <casacore/measures/Measures/MPosition.h>
92 : #include <casacore/casa/Quanta/MVEpoch.h>
93 : #include <casacore/casa/Quanta/MVTime.h>
94 : #include <casacore/measures/Measures/MEpoch.h>
95 : #include <casacore/measures/Measures/MeasTable.h>
96 : #include <casacore/measures/Measures/MeasFrame.h>
97 :
98 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
99 : #include <casacore/ms/MeasurementSets/MSColumns.h>
100 : #include <casacore/ms/MSSel/MSSelection.h>
101 : #include <casacore/ms/MSSel/MSDataDescIndex.h>
102 : #include <casacore/ms/MeasurementSets/MSDopplerUtil.h>
103 : #include <casacore/ms/MSSel/MSSourceIndex.h>
104 : #include <casacore/ms/MSOper/MSSummary.h>
105 : #include <synthesis/MeasurementEquations/CubeSkyEquation.h>
106 : // Disabling Imager::correct() (gmoellen 06Nov20)
107 : //#include <synthesis/MeasurementEquations/VisEquation.h>
108 : #include <synthesis/MeasurementComponents/ImageSkyModel.h>
109 : #include <synthesis/MeasurementComponents/CEMemImageSkyModel.h>
110 : #include <synthesis/MeasurementComponents/MFCEMemImageSkyModel.h>
111 : #include <synthesis/MeasurementComponents/MFCleanImageSkyModel.h>
112 : #include <synthesis/MeasurementComponents/CSCleanImageSkyModel.h>
113 : #include <synthesis/MeasurementComponents/MFMSCleanImageSkyModel.h>
114 : #include <synthesis/MeasurementComponents/HogbomCleanImageSkyModel.h>
115 : #include <synthesis/MeasurementComponents/MSCleanImageSkyModel.h>
116 : #include <synthesis/MeasurementComponents/NNLSImageSkyModel.h>
117 : #include <synthesis/MeasurementComponents/WBCleanImageSkyModel.h>
118 : #include <synthesis/MeasurementComponents/MultiThreadedVisResampler.h>
119 : #include <synthesis/MeasurementComponents/GridBoth.h>
120 : #include <synthesis/TransformMachines/rGridFT.h>
121 : #include <synthesis/TransformMachines/SetJyGridFT.h>
122 : #include <synthesis/TransformMachines/MosaicFT.h>
123 : #include <synthesis/TransformMachines/WProjectFT.h>
124 : #include <synthesis/MeasurementComponents/nPBWProjectFT.h>
125 : #include <synthesis/MeasurementComponents/PBMosaicFT.h>
126 : #include <synthesis/TransformMachines/PBMath.h>
127 : #include <synthesis/TransformMachines/SimpleComponentFTMachine.h>
128 : #include <synthesis/TransformMachines/SimpCompGridMachine.h>
129 : #include <synthesis/TransformMachines/VPSkyJones.h>
130 : #include <synthesis/TransformMachines/SynthesisError.h>
131 : #include <synthesis/TransformMachines/HetArrayConvFunc.h>
132 : #include <synthesis/TransformMachines/VisibilityResamplerBase.h>
133 :
134 : #include <synthesis/DataSampling/SynDataSampling.h>
135 : #include <synthesis/DataSampling/SDDataSampling.h>
136 : #include <synthesis/DataSampling/ImageDataSampling.h>
137 :
138 : #include <synthesis/TransformMachines/StokesImageUtil.h>
139 : #include <casacore/lattices/LRegions/LattRegionHolder.h>
140 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
141 : #include <casacore/lattices/Lattices/LatticeIterator.h>
142 : #include <casacore/lattices/LEL/LatticeExpr.h>
143 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
144 : #include <casacore/lattices/LRegions/LCEllipsoid.h>
145 : #include <casacore/lattices/LRegions/LCRegion.h>
146 : #include <casacore/lattices/LRegions/LCBox.h>
147 : #include <casacore/lattices/LRegions/LCIntersection.h>
148 : #include <casacore/lattices/LRegions/LCUnion.h>
149 : #include <casacore/lattices/LRegions/LCExtension.h>
150 :
151 : #include <casacore/images/Images/ImageRegrid.h>
152 : #include <casacore/images/Regions/ImageRegion.h>
153 : #include <casacore/images/Regions/RegionManager.h>
154 : #include <casacore/images/Regions/WCBox.h>
155 : #include <casacore/images/Regions/WCUnion.h>
156 : #include <casacore/images/Regions/WCIntersection.h>
157 : #include <synthesis/TransformMachines/PBMath.h>
158 : #include <casacore/images/Images/PagedImage.h>
159 : #include <casacore/images/Images/ImageInfo.h>
160 : #include <casacore/images/Images/SubImage.h>
161 : #include <casacore/images/Images/ImageUtilities.h>
162 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
163 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
164 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
165 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
166 : #include <casacore/coordinates/Coordinates/Projection.h>
167 : #include <casacore/coordinates/Coordinates/ObsInfo.h>
168 :
169 : #include <components/ComponentModels/ComponentList.h>
170 : #include <components/ComponentModels/ConstantSpectrum.h>
171 : #include <components/ComponentModels/TabularSpectrum.h>
172 : #include <components/ComponentModels/Flux.h>
173 : #include <components/ComponentModels/FluxStandard.h>
174 : #include <components/ComponentModels/PointShape.h>
175 : #include <components/ComponentModels/DiskShape.h>
176 :
177 : #include <casacore/casa/OS/HostInfo.h>
178 :
179 : #include <components/ComponentModels/ComponentList.h>
180 :
181 : #include <casacore/measures/Measures/UVWMachine.h>
182 :
183 : #include <sstream>
184 :
185 : #include <sys/types.h>
186 : #include <unistd.h>
187 : #include <limits>
188 :
189 : #include <synthesis/TransformMachines/AWProjectFT.h>
190 : #include <synthesis/TransformMachines/AWProjectWBFT.h>
191 : #include <synthesis/TransformMachines/MultiTermFT.h>
192 : #include <synthesis/TransformMachines/NewMultiTermFT.h>
193 : #include <synthesis/TransformMachines/AWConvFunc.h>
194 : //#include <synthesis/TransformMachines/AWConvFuncEPJones.h>
195 : #include <synthesis/TransformMachines/NoOpATerm.h>
196 :
197 : #include <synthesis/Utilities/PointingDirectionCalculator.h>
198 : #include <synthesis/Utilities/SingleDishBeamUtil.h>
199 :
200 : using namespace std;
201 :
202 : using namespace casacore;
203 : namespace casa { //# NAMESPACE CASA - BEGIN
204 :
205 0 : String Imager::imageName()
206 : {
207 0 : LogIO os(LogOrigin("imager", "imageName()", WHERE));
208 : try {
209 0 : lock();
210 0 : String name(msname_p);
211 0 : MSColumns msc(*ms_p);
212 0 : if(datafieldids_p.shape() !=0) {
213 0 : name=msc.field().name()(datafieldids_p(0));
214 : }
215 0 : else if(fieldid_p > -1) {
216 0 : name=msc.field().name()(fieldid_p);
217 : }
218 0 : unlock();
219 0 : return name;
220 0 : } catch (AipsError x) {
221 0 : unlock();
222 0 : os << LogIO::SEVERE << "Caught Exception: "<< x.getMesg() << LogIO::EXCEPTION;
223 0 : return "";
224 0 : }
225 : return String("imagerImage");
226 0 : }
227 :
228 : // imagecoordinates2 (use subMS method to get freq vectors)
229 : // Make standard choices for coordinates
230 0 : Bool Imager::imagecoordinates2(CoordinateSystem& coordInfo, const Bool verbose)
231 : {
232 0 : if(!valid()) return false;
233 0 : if(!assertDefinedImageParameters()) return false;
234 0 : LogIO os(LogOrigin("Imager", "imagecoordinates()", WHERE));
235 :
236 : //===Adjust for some setimage defaults if imageNchan=-1 and spwids=-1
237 0 : if(imageNchan_p <=0){
238 0 : imageNchan_p=1;
239 : }
240 0 : if(spectralwindowids_p.nelements()==1){
241 0 : if(spectralwindowids_p[0]<0){
242 0 : spectralwindowids_p.resize();
243 0 : if(dataspectralwindowids_p.nelements()==0){
244 0 : Int nspwinms=ms_p->spectralWindow().nrow();
245 0 : dataspectralwindowids_p.resize(nspwinms);
246 0 : indgen(dataspectralwindowids_p);
247 : }
248 0 : spectralwindowids_p=dataspectralwindowids_p;
249 : }
250 :
251 : }
252 0 : if(fieldid_p < 0){
253 0 : if(datafieldids_p.shape() !=0) {
254 0 : fieldid_p=datafieldids_p[0];
255 : }
256 : else{
257 0 : fieldid_p=0; //best default if nothing is specified
258 : }
259 : }
260 : //===end of default
261 0 : Vector<Double> deltas(2);
262 0 : deltas(0)=-mcellx_p.get("rad").getValue();
263 0 : deltas(1)=mcelly_p.get("rad").getValue();
264 :
265 : //MSColumns msc(*ms_p);
266 0 : MSColumns msc(*mssel_p); // CAS-11503
267 0 : MFrequency::Types obsFreqRef=MFrequency::DEFAULT;
268 0 : ScalarColumn<Int> measFreqRef(ms_p->spectralWindow(),
269 0 : MSSpectralWindow::columnName(MSSpectralWindow::MEAS_FREQ_REF));
270 : //using the first frame of reference; TO DO should do the right thing
271 : //for different frames selected.
272 : //Int eh = spectralwindowids_p(0);
273 0 : if(spectralwindowids_p.size() && measFreqRef(spectralwindowids_p(0)) >=0) {
274 0 : obsFreqRef=(MFrequency::Types)measFreqRef(spectralwindowids_p(0));
275 : }
276 :
277 :
278 0 : MVDirection mvPhaseCenter(phaseCenter_p.getAngle());
279 : // Normalize correctly
280 0 : MVAngle ra=mvPhaseCenter.get()(0);
281 0 : ra(0.0);
282 0 : MVAngle dec=mvPhaseCenter.get()(1);
283 0 : Vector<Double> refCoord(2);
284 0 : refCoord(0)=ra.get().getValue();
285 0 : refCoord(1)=dec;
286 :
287 0 : Vector<Double> refPixel(2);
288 0 : refPixel(0) = Double(nx_p / 2);
289 0 : refPixel(1) = Double(ny_p / 2);
290 :
291 : //defining observatory...needed for position on earth
292 0 : String telescop = msc.observation().telescopeName()(0);
293 :
294 : // defining epoch as begining time from timerange in OBSERVATION subtable
295 : // Using first observation for now
296 : //MEpoch obsEpoch = msc.observation().timeRangeMeas()(0)(IPosition(1,0));
297 : // modified to use main table's TIME column for better match with what
298 : // VisIter does.
299 0 : MEpoch obsEpoch = msc.timeMeas()(0);
300 :
301 : //Now finding the position of the telescope on Earth...needed for proper
302 : //frequency conversions
303 :
304 0 : MPosition obsPosition;
305 0 : if(nonDefaultLocation()) {
306 : os << LogIO::DEBUG1 << "Using user defined location: "
307 0 : << mLocation_p.getRefString() << " " << mLocation_p.get("m")
308 0 : << LogIO::POST;
309 0 : obsPosition = mLocation_p;
310 0 : freqFrameValid_p = true;
311 0 : } else if(!(MeasTable::Observatory(obsPosition, telescop))){
312 : os << LogIO::WARN << "Did not get the position of " << telescop
313 0 : << " from data repository" << LogIO::POST;
314 : os << LogIO::WARN
315 : << "Please contact CASA to add it to the repository."
316 0 : << LogIO::POST;
317 0 : os << LogIO::WARN << "Frequency conversion will not work " << LogIO::POST;
318 0 : freqFrameValid_p = false;
319 : }
320 : else{
321 0 : mLocation_p = obsPosition;
322 0 : freqFrameValid_p = true;
323 : os << LogIO::DEBUG1 << "Using observatory location of "
324 0 : << telescop << ": " << mLocation_p.getRefString()
325 0 : << " " << mLocation_p.get("m") << LogIO::POST;
326 : }
327 : //Make sure frame conversion is switched off for REST frame data.
328 0 : freqFrameValid_p=freqFrameValid_p && (obsFreqRef !=MFrequency::REST);
329 :
330 : // Now find the projection to use: could probably also use
331 : // max(abs(w))=0.0 as a criterion
332 0 : Projection::Type ptype = Projection::type(projection_p);
333 0 : Projection projection(ptype);
334 0 : if(ptype == Projection::SIN
335 0 : && (telescop == "ATCASCP" || telescop == "WSRT" || telescop == "DRAO")) {
336 : os << LogIO::NORMAL // Loglevel NORMAL
337 : << "Using SIN image projection adjusted for "
338 0 : << (telescop == "ATCASCP" ? 'S' : 'N') << "CP"
339 0 : << LogIO::POST;
340 0 : Vector<Double> projectionParameters(2);
341 0 : projectionParameters(0) = 0.0;
342 0 : if(sin(dec) != 0.0){
343 0 : projectionParameters(1) = cos(dec)/sin(dec);
344 0 : projection=Projection(Projection::SIN, projectionParameters);
345 : }
346 : else {
347 : os << LogIO::WARN
348 : << "Singular projection for " << telescop << ": using plain SIN"
349 0 : << LogIO::POST;
350 0 : projection=Projection(Projection::SIN);
351 : }
352 0 : }
353 : else {
354 0 : os << LogIO::DEBUGGING << "Using " << projection_p << " image projection" << LogIO::POST;
355 : }
356 0 : os << LogIO::NORMAL;
357 :
358 0 : Matrix<Double> xform(2,2);
359 0 : xform=0.0;xform.diagonal()=1.0;
360 : DirectionCoordinate
361 0 : myRaDec(MDirection::Types(phaseCenter_p.getRefPtr()->getType()),
362 : projection,
363 0 : refCoord(0), refCoord(1),
364 0 : deltas(0), deltas(1),
365 : xform,
366 0 : refPixel(0), refPixel(1));
367 :
368 : // Now set up spectral coordinate
369 0 : SpectralCoordinate* mySpectral=0;
370 0 : Double refChan=0.0;
371 :
372 : // Spectral synthesis
373 : // For mfs band we set the window to include all spectral windows
374 0 : Int nspw=spectralwindowids_p.nelements();
375 0 : if (imageMode_p=="MFS") {
376 0 : Double fmin=C::dbl_max;
377 0 : Double fmax=-(C::dbl_max);
378 0 : Double fmean=0.0;
379 : /*
380 : Int nms = freqrange_p.shape()(0);
381 : for (Int i=0;i<nspw;++i) {
382 : Int spw=spectralwindowids_p(i);
383 : Vector<Double> chanFreq=msc.spectralWindow().chanFreq()(spw);
384 : Vector<Double> freqResolution = msc.spectralWindow().chanWidth()(spw);
385 :
386 : if(dataMode_p=="none"){
387 :
388 : if(nms >1) {
389 : for (Int j=0; j<nms; ++j) {
390 : fmin=min(fmin,freqrange_p(j,0));
391 : fmax=max(fmax,freqrange_p(j,1));
392 : }
393 : }
394 : else if(i==0) {
395 : fmin=min(chanFreq-abs(0.5*freqResolution));
396 : fmax=max(chanFreq+abs(0.5*freqResolution));
397 : }
398 : else {
399 : fmin=min(fmin,min(chanFreq-abs(0.5*freqResolution)));
400 : fmax=max(fmax,max(chanFreq+abs(0.5*freqResolution)));
401 : }
402 : }
403 : else if(dataMode_p=="channel"){
404 : // This needs some careful thought about roundoff - it is likely
405 : // still adding an extra half-channel at top and bottom but
406 : // if the freqResolution is nonlinear, there are subtleties
407 : // multiple MS case
408 : if(nms >1) {
409 : for (Int j=0; j<nms; ++j) {
410 : fmin=min(fmin,freqrange_p(j,0));
411 : fmax=max(fmax,freqrange_p(j,1));
412 : }
413 : }
414 : // single ms case
415 : else {
416 : Int elnchan=chanFreq.nelements();
417 : Int firstchan=0;
418 : Int elstep=1;
419 : for (uInt jj=0; jj < dataspectralwindowids_p.nelements(); ++jj){
420 : if(dataspectralwindowids_p[jj]==spw){
421 : firstchan=dataStart_p[jj];
422 : elnchan=dataNchan_p[jj];
423 : elstep=dataStep_p[jj];
424 : }
425 : }
426 : Int lastchan=firstchan+ elnchan*elstep;
427 : for(Int k=firstchan ; k < lastchan ; k+=elstep){
428 : fmin=min(fmin,chanFreq[k]-abs(freqResolution[k]*(elstep-0.5)));
429 : fmax=max(fmax,chanFreq[k]+abs(freqResolution[k]*(elstep-0.5)));
430 : }
431 : }
432 : }
433 : else{
434 : this->unlock();
435 : os << LogIO::SEVERE
436 : << "setdata has to be in 'channel' or 'none' mode for 'mfs' imaging to work"
437 : << LogIO::EXCEPTION;
438 : return false;
439 : }
440 :
441 : }
442 :
443 : */
444 0 : rvi_p->getFreqInSpwRange(fmin, fmax, freqFrameValid_p ? MFrequency::LSRK : freqFrame_p);
445 :
446 :
447 0 : fmean=(fmax+fmin)/2.0;
448 0 : Vector<Double> restFreqArray;
449 0 : Double restFreq=fmean;
450 0 : if(getRestFreq(restFreqArray, spectralwindowids_p(0))){
451 0 : restFreq=restFreqArray[0];
452 : }
453 0 : imageNchan_p=1;
454 0 : Double finc=(fmax-fmin);
455 0 : mySpectral = new SpectralCoordinate(freqFrameValid_p ? MFrequency::LSRK : freqFrame_p, fmean//-finc/2.0
456 : , finc,
457 0 : refChan, restFreq);
458 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3) // Loglevel INFO
459 : << "Center frequency = "
460 0 : << MFrequency(Quantity(fmean, "Hz")).get("GHz").getValue()
461 : << " GHz, synthesized continuum bandwidth = "
462 0 : << MFrequency(Quantity(finc, "Hz")).get("GHz").getValue()
463 0 : << " GHz" << LogIO::POST;
464 :
465 0 : if(ntaylor_p>1 && reffreq_p==0.0)
466 : {
467 0 : reffreq_p = fmean;
468 0 : os << "Setting center frequency as MS-MFS reference frequency" << LogIO::POST;
469 : }
470 0 : }
471 :
472 0 : else if(imageMode_p.contains("FREQ")) {
473 0 : if(imageNchan_p==0) {
474 0 : this->unlock();
475 : os << LogIO::SEVERE << "Must specify number of channels"
476 0 : << LogIO::EXCEPTION;
477 0 : return false;
478 : }
479 0 : Double restFreq=mfImageStart_p.get("Hz").getValue();
480 0 : Vector<Double> restFreqVec;
481 0 : if(getRestFreq(restFreqVec, spectralwindowids_p(0))){
482 0 : restFreq=restFreqVec[0];
483 : }
484 0 : MFrequency::Types mfreqref=(obsFreqRef==(MFrequency::REST)) ? MFrequency::REST : MFrequency::castType(mfImageStart_p.getRef().getType()) ;
485 :
486 : /////Some problem here it is really goofing up in getting frequency
487 : // -> fixed was made in calcImFreqs -TT
488 0 : Vector<Double> imgridfreqs;
489 0 : Vector<Double> imfreqres;
490 : //rstate=calcImFreqs(imgridfreqs, imfreqres, mfreqref, obsEpoch, obsPosition,restFreq);
491 : // should use obsFreqRef
492 :
493 0 : if (!calcImFreqs(imgridfreqs, imfreqres, obsFreqRef, obsEpoch, mLocation_p ,restFreq)) {
494 0 : return false;
495 : }
496 : //cerr<<"imfreqres(0)="<<imfreqres(0)<<endl;
497 :
498 :
499 0 : if (imageNchan_p==1) {
500 0 : mySpectral = new SpectralCoordinate(mfreqref,
501 0 : mfImageStart_p.get("Hz").getValue(),
502 0 : mfImageStep_p.get("Hz").getValue(),
503 0 : refChan, restFreq);
504 : }
505 : else {
506 0 : Double finc= imgridfreqs(1)-imgridfreqs(0);
507 0 : mySpectral = new SpectralCoordinate(mfreqref, imgridfreqs(0), finc, refChan, restFreq);
508 : //cerr << "after myspectral2 " << mySpectral->referenceValue() << " pixel " << mySpectral->referencePixel() << endl;
509 : //debug TT
510 : //Double wrld,pixl;
511 : //pixl=0.0;
512 : //mySpectral->toWorld(wrld,pixl);
513 : //cerr<<"world="<<wrld<<" pixel="<<pixl;
514 : }
515 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
516 : << "Start frequency = " // Loglevel INFO
517 0 : << mfImageStart_p.get("GHz").getValue()
518 : << ", channel increment = "
519 0 : << mfImageStep_p.get("GHz").getValue()
520 : << "GHz, frequency frame = "
521 : << MFrequency::showType(mfreqref)
522 0 : << endl;
523 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
524 : << "Rest frequency is " // Loglevel INFO
525 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
526 0 : << "GHz" << LogIO::POST;
527 :
528 0 : }
529 :
530 :
531 : else {
532 0 : Vector<Double> chanFreq;
533 0 : Vector<Double> freqResolution;
534 : //starting with a default rest frequency to be ref
535 : //in case none is defined
536 : Double restFreq=
537 0 : msc.spectralWindow().refFrequency()(spectralwindowids_p(0));
538 :
539 0 : for (Int spwIndex=0; spwIndex < nspw; ++spwIndex){
540 :
541 0 : Int spw=spectralwindowids_p(spwIndex);
542 :
543 0 : Vector<Double> restFreqArray;
544 0 : if(getRestFreq(restFreqArray, spw)){
545 0 : if(spwIndex==0){
546 0 : restFreq = restFreqArray(0);
547 : }
548 : else{
549 0 : if(restFreq != restFreqArray(0)){
550 : os << LogIO::WARN << "Rest frequencies are different for spectralwindows selected "
551 0 : << LogIO::POST;
552 : os << LogIO::WARN
553 : <<"Will be using the restFreq defined in spectralwindow "
554 0 : << spectralwindowids_p(0) << LogIO::POST;
555 : }
556 :
557 : }
558 : }
559 0 : }
560 :
561 : // use data frame here (for channel mode)
562 0 : if (!calcImFreqs(chanFreq, freqResolution, obsFreqRef, obsEpoch, mLocation_p,restFreq)) {
563 0 : return false;
564 : }
565 :
566 0 : if(imageMode_p=="CHANNEL") {
567 0 : if(imageNchan_p==0) {
568 0 : this->unlock();
569 : os << LogIO::SEVERE << "Must specify number of channels"
570 0 : << LogIO::EXCEPTION;
571 0 : return false;
572 : }
573 0 : if(imageStep_p==0)
574 0 : imageStep_p=1;
575 : // TT: commented these out otherwise the case for multiple MSes would not work
576 : // Int nsubchans=
577 : // (chanFreq.shape()(0) - Int(imageStart_p)+1)/Int(imageStep_p)+1;
578 : // if((nsubchans >0) && (imageNchan_p>nsubchans)) imageNchan_p=nsubchans;
579 :
580 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
581 : << "Image spectral coordinate: "<< imageNchan_p
582 : << " channels, starting at visibility channel "
583 0 : << imageStart_p << " stepped by " << imageStep_p << LogIO::POST;
584 :
585 : Double finc;
586 : {
587 : // Now use outframe (instead of data frame) as the rest of
588 : // the modes do
589 : //
590 :
591 0 : finc= ((chanFreq.shape().nelements()==1) && (chanFreq.shape()[0] > 1)) ? chanFreq(1)-chanFreq(0): freqResolution[0];
592 :
593 0 : mySpectral = new SpectralCoordinate(freqFrame_p,
594 0 : chanFreq(0),
595 : finc,
596 0 : refChan, restFreq);
597 : }
598 :
599 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
600 : << "Frequency = " // Loglevel INFO
601 0 : << MFrequency(Quantity(chanFreq(0), "Hz")).get("GHz").getValue()
602 : << ", channel increment = "
603 0 : << MFrequency(Quantity(finc, "Hz")).get("GHz").getValue()
604 0 : << "GHz" << endl;
605 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
606 : << "Rest frequency is " // Loglevel INFO
607 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
608 0 : << "GHz" << LogIO::POST;
609 :
610 : }
611 : // Spectral channels resampled at equal increments in optical velocity
612 : // Here we compute just the first two channels and use increments for
613 : // the others
614 0 : else if (imageMode_p=="VELOCITY" || imageMode_p.contains("RADIO")) {
615 0 : if(imageNchan_p==0) {
616 0 : this->unlock();
617 : os << LogIO::SEVERE << "Must specify number of channels"
618 0 : << LogIO::EXCEPTION;
619 0 : return false;
620 : }
621 : {
622 0 : ostringstream oos;
623 0 : oos << "Image spectral coordinate:"<< imageNchan_p
624 0 : << " channels, starting at radio velocity " << mImageStart_p
625 0 : << " stepped by " << mImageStep_p << endl;
626 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
627 0 : << String(oos); // Loglevel INFO
628 0 : }
629 :
630 0 : MFrequency::Types mfreqref=MFrequency::LSRK;
631 : //Can't convert to frame in mImageStart
632 0 : if(!MFrequency::getType(mfreqref, (MRadialVelocity::showType(mImageStart_p.getRef().getType()))))
633 0 : mfreqref=freqFrame_p;
634 0 : mfreqref=(obsFreqRef==(MFrequency::REST)) ? MFrequency::REST : mfreqref;
635 : //rstate=calcImFreqs(chanFreq, freqResolution, mfreqref, obsEpoch, obsPosition,restFreq);
636 0 : if (!calcImFreqs(chanFreq, freqResolution, obsFreqRef, obsEpoch, mLocation_p, restFreq)) {
637 0 : return false;
638 : }
639 :
640 :
641 : Double finc;
642 0 : if(imageNchan_p ==1) {
643 0 : finc = freqResolution(0);
644 : }
645 : else {
646 0 : finc = chanFreq(1)-chanFreq(0);
647 : }
648 : //mySpectral = new SpectralCoordinate(obsFreqRef,
649 0 : mySpectral = new SpectralCoordinate(mfreqref,
650 0 : chanFreq(0),
651 : finc,
652 0 : refChan, restFreq);
653 :
654 : {
655 0 : ostringstream oos;
656 0 : oos << "Reference Frequency = "
657 0 : << MFrequency(Quantity(chanFreq(0), "Hz")).get("GHz")
658 0 : << ", spectral increment = "
659 0 : << MFrequency(Quantity(finc, "Hz")).get("GHz")
660 0 : << ", frequency frame = "
661 0 : << MFrequency::showType(mfreqref)
662 0 : << endl;
663 0 : oos << "Rest frequency is "
664 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
665 0 : << " GHz" << endl;
666 0 : os << LogIO::NORMAL << String(oos) << LogIO::POST; // Loglevel INFO
667 0 : }
668 :
669 : }
670 : // Since optical/relativistic velocity is non-linear in frequency, we have to
671 : // pass in all the frequencies. For radio velocity we can use
672 : // a linear axis.
673 0 : else if (imageMode_p=="OPTICALVELOCITY" || imageMode_p.contains("OPTICAL") || imageMode_p.contains("TRUE")
674 0 : || imageMode_p.contains("BETA") || imageMode_p.contains("RELATI") ) {
675 0 : if(imageNchan_p==0) {
676 0 : this->unlock();
677 : os << LogIO::SEVERE << "Must specify number of channels"
678 0 : << LogIO::EXCEPTION;
679 0 : return false;
680 : }
681 : {
682 0 : ostringstream oos;
683 0 : oos << "Image spectral coordinate: "<< imageNchan_p
684 0 : << " channels, starting at optical velocity " << mImageStart_p
685 0 : << " stepped by " << mImageStep_p << endl;
686 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
687 0 : << String(oos); // Loglevel INFO
688 0 : }
689 :
690 : // Use this next line when non-linear is working
691 : // when selecting in velocity its specfied freqframe or REST
692 0 : MFrequency::Types imfreqref=(obsFreqRef==MFrequency::REST) ? MFrequency::REST : freqFrame_p;
693 : //rstate=calcImFreqs(chanFreq, freqResolution, imfreqref, obsEpoch, obsPosition,restFreq);
694 0 : if (!calcImFreqs(chanFreq, freqResolution, obsFreqRef, obsEpoch, mLocation_p,restFreq)) {
695 0 : return false;
696 : }
697 :
698 0 : if (imageNchan_p==1) {
699 0 : mySpectral = new SpectralCoordinate(imfreqref,
700 0 : chanFreq(0),
701 0 : freqResolution(0),
702 0 : refChan, restFreq);
703 : }
704 : else {
705 0 : mySpectral = new SpectralCoordinate(imfreqref, chanFreq, restFreq);
706 : }
707 :
708 : {
709 0 : ostringstream oos;
710 0 : oos << "Reference Frequency = "
711 : //<< MFrequency(Quantity(freqs(0), "Hz")).get("GHz")
712 0 : << MFrequency(Quantity(chanFreq(0), "Hz")).get("GHz")
713 : << " Ghz, "
714 0 : <<" frequency frame= "<<MFrequency::showType(imfreqref)<<endl;
715 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
716 0 : << String(oos) << LogIO::POST; // Loglevel INFO
717 0 : }
718 : }
719 : else {
720 0 : this->unlock();
721 0 : os << LogIO::SEVERE << "Unknown mode " << imageMode_p
722 0 : << LogIO::EXCEPTION;
723 0 : return false;
724 : }
725 :
726 :
727 0 : }
728 :
729 : //In FTMachine lsrk is used for channel matching with data channel
730 : //hence we make sure that
731 : // we convert to lsrk when dealing with the channels
732 0 : freqFrameValid_p=freqFrameValid_p && (obsFreqRef !=MFrequency::REST);
733 0 : if(freqFrameValid_p){
734 0 : mySpectral->setReferenceConversion(MFrequency::LSRK, obsEpoch,
735 : obsPosition,
736 0 : phaseCenter_p);
737 : }
738 :
739 : // Polarization
740 0 : Vector<String> polType=msc.feed().polarizationType()(0);
741 0 : if (polType(0)!="X" && polType(0)!="Y" &&
742 0 : polType(0)!="R" && polType(0)!="L") {
743 : os << LogIO::WARN << "Unknown stokes types in feed table: ["
744 0 : << polType(0) << ", " << polType(1) << "]" << endl
745 0 : << "Results open to question!" << LogIO::POST;
746 : }
747 :
748 0 : if (polType(0)=="X" || polType(0)=="Y") {
749 0 : polRep_p=StokesImageUtil::LINEAR;
750 : os << LogIO::DEBUG1
751 0 : << "Preferred polarization representation is linear" << LogIO::POST;
752 : }
753 : else {
754 0 : polRep_p=StokesImageUtil::CIRCULAR;
755 : os << LogIO::DEBUG1
756 0 : << "Preferred polarization representation is circular" << LogIO::POST;
757 : }
758 :
759 : // Compare user input with whatever is allowed by the data.
760 : // If possible, allow.
761 :
762 0 : Vector<Int> whichStokes = decideNPolPlanes(true);
763 0 : if(whichStokes.nelements()==0 || (whichStokes.nelements()==1 && whichStokes[0]==0) )
764 : {
765 0 : if(polRep_p==StokesImageUtil::CIRCULAR)
766 0 : os << LogIO::SEVERE << "Stokes selection of " << stokes_p << " is not valid for Circular feeds." << LogIO::EXCEPTION;
767 : else
768 0 : os << LogIO::SEVERE << "Stokes selection of " << stokes_p << " is not valid for Linear feeds." << LogIO::EXCEPTION;
769 0 : return false;
770 : }
771 :
772 0 : StokesCoordinate myStokes(whichStokes);
773 :
774 : // os << LogIO::DEBUG1 << "imagecoordinate : " << (coordInfo).stokesCoordinate((coordInfo).findCoordinate(Coordinate::STOKES)).stokes() << LogIO::POST;
775 :
776 : //Set Observatory info
777 0 : ObsInfo myobsinfo;
778 0 : myobsinfo.setTelescope(telescop);
779 0 : myobsinfo.setPointingCenter(mvPhaseCenter);
780 0 : myobsinfo.setObsDate(obsEpoch);
781 0 : myobsinfo.setObserver(msc.observation().observer()(0));
782 0 : this->setObsInfo(myobsinfo);
783 :
784 : //Adding everything to the coordsystem
785 0 : coordInfo.addCoordinate(myRaDec);
786 0 : coordInfo.addCoordinate(myStokes);
787 0 : coordInfo.addCoordinate(*mySpectral);
788 0 : coordInfo.setObsInfo(myobsinfo);
789 :
790 0 : if(mySpectral) delete mySpectral;
791 :
792 0 : return true;
793 0 : }
794 :
795 : //Return if mLocation_p is non-default MPosition
796 0 : Bool Imager::nonDefaultLocation(){
797 0 : MPosition default_MPosition;
798 : //Check for measure reference
799 0 : if (mLocation_p.getRefString() != default_MPosition.getRefString())
800 0 : return true;
801 : //Check for measure value
802 0 : if (mLocation_p.get("m") != default_MPosition.get("m"))
803 0 : return true;
804 0 : return false;
805 0 : }
806 :
807 :
808 : // Make standard choices for coordinates
809 0 : Bool Imager::imagecoordinates(CoordinateSystem& coordInfo, const Bool verbose)
810 : {
811 0 : if(!valid()) return false;
812 0 : if(!assertDefinedImageParameters()) return false;
813 0 : LogIO os(LogOrigin("Imager", "imagecoordinates()", WHERE));
814 :
815 : //===Adjust for some setimage defaults if imageNchan=-1 and spwids=-1
816 0 : if(imageNchan_p <=0){
817 0 : imageNchan_p=1;
818 : }
819 0 : if(spectralwindowids_p.nelements()==1){
820 0 : if(spectralwindowids_p[0]<0){
821 0 : spectralwindowids_p.resize();
822 0 : if(dataspectralwindowids_p.nelements()==0){
823 0 : Int nspwinms=ms_p->spectralWindow().nrow();
824 0 : dataspectralwindowids_p.resize(nspwinms);
825 0 : indgen(dataspectralwindowids_p);
826 : }
827 0 : spectralwindowids_p=dataspectralwindowids_p;
828 : }
829 :
830 : }
831 0 : if(fieldid_p < 0){
832 0 : if(datafieldids_p.shape() !=0) {
833 0 : fieldid_p=datafieldids_p[0];
834 : }
835 : else{
836 0 : fieldid_p=0; //best default if nothing is specified
837 : }
838 : }
839 : //===end of default
840 0 : Vector<Double> deltas(2);
841 0 : deltas(0)=-mcellx_p.get("rad").getValue();
842 0 : deltas(1)=mcelly_p.get("rad").getValue();
843 :
844 0 : MSColumns msc(*ms_p);
845 0 : MFrequency::Types obsFreqRef=MFrequency::DEFAULT;
846 0 : ScalarColumn<Int> measFreqRef(ms_p->spectralWindow(),
847 0 : MSSpectralWindow::columnName(MSSpectralWindow::MEAS_FREQ_REF));
848 : //using the first frame of reference; TO DO should do the right thing
849 : //for different frames selected.
850 : //Int eh = spectralwindowids_p(0);
851 0 : if(spectralwindowids_p.size() && measFreqRef(spectralwindowids_p(0)) >=0) {
852 0 : obsFreqRef=(MFrequency::Types)measFreqRef(spectralwindowids_p(0));
853 : }
854 :
855 :
856 0 : MVDirection mvPhaseCenter(phaseCenter_p.getAngle());
857 : // Normalize correctly
858 0 : MVAngle ra=mvPhaseCenter.get()(0);
859 0 : ra(0.0);
860 0 : MVAngle dec=mvPhaseCenter.get()(1);
861 0 : Vector<Double> refCoord(2);
862 0 : refCoord(0)=ra.get().getValue();
863 0 : refCoord(1)=dec;
864 :
865 0 : Vector<Double> refPixel(2);
866 0 : refPixel(0) = Double(nx_p / 2);
867 0 : refPixel(1) = Double(ny_p / 2);
868 :
869 : //defining observatory...needed for position on earth
870 0 : String telescop = msc.observation().telescopeName()(0);
871 :
872 : // defining epoch as begining time from timerange in OBSERVATION subtable
873 : // Using first observation for now
874 : //MEpoch obsEpoch = msc.observation().timeRangeMeas()(0)(IPosition(1,0));
875 : // modified to use main table's TIME column for better match with what
876 : // VisIter does.
877 0 : MEpoch obsEpoch = msc.timeMeas()(0);
878 :
879 : //Now finding the position of the telescope on Earth...needed for proper
880 : //frequency conversions
881 :
882 0 : MPosition obsPosition;
883 0 : if(nonDefaultLocation()) {
884 0 : freqFrameValid_p = true;
885 0 : } else if(!(MeasTable::Observatory(obsPosition, telescop))){
886 : os << LogIO::WARN << "Did not get the position of " << telescop
887 0 : << " from data repository" << LogIO::POST;
888 : os << LogIO::WARN
889 : << "Please contact CASA to add it to the repository."
890 0 : << LogIO::POST;
891 0 : os << LogIO::WARN << "Frequency conversion will not work " << LogIO::POST;
892 0 : freqFrameValid_p = false;
893 : }
894 : else{
895 0 : mLocation_p = obsPosition;
896 0 : freqFrameValid_p = true;
897 : }
898 : //Make sure frame conversion is switched off for REST frame data.
899 0 : freqFrameValid_p=freqFrameValid_p && (obsFreqRef !=MFrequency::REST);
900 : // Now find the projection to use: could probably also use
901 : // max(abs(w))=0.0 as a criterion
902 0 : Projection projection(Projection::SIN);
903 0 : if(telescop == "ATCASCP" || telescop == "WSRT" || telescop == "DRAO") {
904 : os << LogIO::NORMAL // Loglevel NORMAL
905 : << "Using SIN image projection adjusted for "
906 0 : << (telescop == "ATCASCP" ? 'S' : 'N') << "CP"
907 0 : << LogIO::POST;
908 0 : Vector<Double> projectionParameters(2);
909 0 : projectionParameters(0) = 0.0;
910 0 : if(sin(dec) != 0.0){
911 0 : projectionParameters(1) = cos(dec)/sin(dec);
912 0 : projection=Projection(Projection::SIN, projectionParameters);
913 : }
914 : else {
915 : os << LogIO::WARN
916 : << "Singular projection for " << telescop << ": using plain SIN"
917 0 : << LogIO::POST;
918 0 : projection=Projection(Projection::SIN);
919 : }
920 0 : }
921 : else {
922 0 : os << LogIO::DEBUGGING << "Using SIN image projection" << LogIO::POST;
923 : }
924 0 : os << LogIO::NORMAL;
925 :
926 0 : Matrix<Double> xform(2,2);
927 0 : xform=0.0;xform.diagonal()=1.0;
928 : DirectionCoordinate
929 0 : myRaDec(MDirection::Types(phaseCenter_p.getRefPtr()->getType()),
930 : projection,
931 0 : refCoord(0), refCoord(1),
932 0 : deltas(0), deltas(1),
933 : xform,
934 0 : refPixel(0), refPixel(1));
935 :
936 : // Now set up spectral coordinate
937 0 : SpectralCoordinate* mySpectral=0;
938 0 : Double refChan=0.0;
939 :
940 : // Spectral synthesis
941 : // For mfs band we set the window to include all spectral windows
942 0 : Int nspw=spectralwindowids_p.nelements();
943 0 : if (imageMode_p=="MFS") {
944 0 : Double fmin=C::dbl_max;
945 0 : Double fmax=-(C::dbl_max);
946 0 : Double fmean=0.0;
947 0 : for (Int i=0;i<nspw;++i) {
948 0 : Int spw=spectralwindowids_p(i);
949 0 : Vector<Double> chanFreq=msc.spectralWindow().chanFreq()(spw);
950 0 : Vector<Double> freqResolution = msc.spectralWindow().chanWidth()(spw);
951 :
952 0 : if(dataMode_p=="none"){
953 :
954 0 : if(i==0) {
955 0 : fmin=min(chanFreq-abs(0.5*freqResolution));
956 0 : fmax=max(chanFreq+abs(0.5*freqResolution));
957 : }
958 : else {
959 0 : fmin=min(fmin,min(chanFreq-abs(0.5*freqResolution)));
960 0 : fmax=max(fmax,max(chanFreq+abs(0.5*freqResolution)));
961 : }
962 : }
963 0 : else if(dataMode_p=="channel"){
964 : // This needs some careful thought about roundoff - it is likely
965 : // still adding an extra half-channel at top and bottom but
966 : // if the freqResolution is nonlinear, there are subtleties
967 0 : Int elnchan=chanFreq.nelements();
968 0 : Int firstchan=0;
969 0 : Int elstep=1;
970 0 : for (uInt jj=0; jj < dataspectralwindowids_p.nelements(); ++jj){
971 0 : if(dataspectralwindowids_p[jj]==spw){
972 0 : firstchan=dataStart_p[jj];
973 0 : elnchan=dataNchan_p[jj];
974 0 : elstep=dataStep_p[jj];
975 : }
976 : }
977 0 : Int lastchan=firstchan+ elnchan*elstep;
978 0 : for(Int k=firstchan ; k < lastchan ; k+=elstep){
979 0 : fmin=min(fmin,chanFreq[k]-abs(freqResolution[k]*(elstep-0.5)));
980 0 : fmax=max(fmax,chanFreq[k]+abs(freqResolution[k]*(elstep-0.5)));
981 : }
982 : }
983 : else{
984 0 : this->unlock();
985 : os << LogIO::SEVERE
986 : << "setdata has to be in 'channel' or 'none' mode for 'mfs' imaging to work"
987 0 : << LogIO::EXCEPTION;
988 0 : return false;
989 : }
990 :
991 0 : }
992 :
993 0 : fmean=(fmax+fmin)/2.0;
994 0 : Vector<Double> restFreqArray;
995 0 : Double restFreq=fmean;
996 0 : if(getRestFreq(restFreqArray, spectralwindowids_p(0))){
997 0 : restFreq=restFreqArray[0];
998 : }
999 0 : imageNchan_p=1;
1000 0 : Double finc=(fmax-fmin);
1001 0 : mySpectral = new SpectralCoordinate(freqFrame_p, fmean//-finc/2.0
1002 : , finc,
1003 0 : refChan, restFreq);
1004 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3) // Loglevel INFO
1005 : << "Center frequency = "
1006 0 : << MFrequency(Quantity(fmean, "Hz")).get("GHz").getValue()
1007 : << " GHz, synthesized continuum bandwidth = "
1008 0 : << MFrequency(Quantity(finc, "Hz")).get("GHz").getValue()
1009 0 : << " GHz" << LogIO::POST;
1010 :
1011 0 : if(ntaylor_p>1 && reffreq_p==0.0)
1012 : {
1013 0 : reffreq_p = fmean;
1014 0 : os << "Setting center frequency as MS-MFS reference frequency" << LogIO::POST;
1015 : }
1016 0 : }
1017 :
1018 0 : else if(imageMode_p.contains("FREQ")) {
1019 0 : if(imageNchan_p==0) {
1020 0 : this->unlock();
1021 : os << LogIO::SEVERE << "Must specify number of channels"
1022 0 : << LogIO::EXCEPTION;
1023 0 : return false;
1024 : }
1025 0 : Double restFreq=mfImageStart_p.get("Hz").getValue();
1026 0 : Vector<Double> restFreqVec;
1027 0 : if(getRestFreq(restFreqVec, spectralwindowids_p(0))){
1028 0 : restFreq=restFreqVec[0];
1029 : }
1030 0 : MFrequency::Types mfreqref=(obsFreqRef==(MFrequency::REST)) ? MFrequency::REST : MFrequency::castType(mfImageStart_p.getRef().getType()) ;
1031 : // mySpectral = new SpectralCoordinate(mfreqref,
1032 : // mfImageStart_p.get("Hz").getValue()+
1033 : //mfImageStep_p.get("Hz").getValue()/2.0,
1034 : // mfImageStep_p.get("Hz").getValue(),
1035 : // refChan, restFreq);
1036 0 : mySpectral = new SpectralCoordinate(mfreqref,
1037 0 : mfImageStart_p.get("Hz").getValue(),
1038 0 : mfImageStep_p.get("Hz").getValue(),
1039 0 : refChan, restFreq);
1040 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1041 : << "Start frequency = " // Loglevel INFO
1042 0 : << mfImageStart_p.get("GHz").getValue()
1043 : << ", channel increment = "
1044 0 : << mfImageStep_p.get("GHz").getValue()
1045 : << "GHz, frequency frame = "
1046 : << MFrequency::showType(mfreqref)
1047 0 : << endl;
1048 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1049 : << "Rest frequency is " // Loglevel INFO
1050 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
1051 0 : << "GHz" << LogIO::POST;
1052 :
1053 0 : }
1054 :
1055 :
1056 : else {
1057 : // if(nspw>1) {
1058 : // os << LogIO::SEVERE << "Line modes allow only one spectral window"
1059 : // << LogIO::POST;
1060 : // return false;
1061 : // }
1062 0 : Vector<Double> chanFreq;
1063 0 : Vector<Double> freqResolution;
1064 : //starting with a default rest frequency to be ref
1065 : //in case none is defined
1066 : Double restFreq=
1067 0 : msc.spectralWindow().refFrequency()(spectralwindowids_p(0));
1068 :
1069 0 : for (Int spwIndex=0; spwIndex < nspw; ++spwIndex){
1070 :
1071 0 : Int spw=spectralwindowids_p(spwIndex);
1072 0 : Int origsize=chanFreq.shape()(0);
1073 0 : Int newsize=origsize+msc.spectralWindow().chanFreq()(spw).shape()(0);
1074 0 : chanFreq.resize(newsize, true);
1075 0 : chanFreq(Slice(origsize, newsize-origsize))=msc.spectralWindow().chanFreq()(spw);
1076 0 : freqResolution.resize(newsize, true);
1077 0 : freqResolution(Slice(origsize, newsize-origsize))=
1078 0 : msc.spectralWindow().chanWidth()(spw);
1079 :
1080 :
1081 :
1082 0 : Vector<Double> restFreqArray;
1083 0 : if(getRestFreq(restFreqArray, spw)){
1084 0 : if(spwIndex==0){
1085 0 : restFreq = restFreqArray(0);
1086 : }
1087 : else{
1088 0 : if(restFreq != restFreqArray(0)){
1089 : os << LogIO::WARN << "Rest frequencies are different for spectralwindows selected "
1090 0 : << LogIO::POST;
1091 : os << LogIO::WARN
1092 : <<"Will be using the restFreq defined in spectralwindow "
1093 0 : << spectralwindowids_p(0) << LogIO::POST;
1094 : }
1095 :
1096 : }
1097 : }
1098 0 : }
1099 :
1100 :
1101 0 : if(imageMode_p=="CHANNEL") {
1102 0 : if(imageNchan_p==0) {
1103 0 : this->unlock();
1104 : os << LogIO::SEVERE << "Must specify number of channels"
1105 0 : << LogIO::EXCEPTION;
1106 0 : return false;
1107 : }
1108 0 : Vector<Double> freqs;
1109 0 : if(imageStep_p==0)
1110 0 : imageStep_p=1;
1111 : // TT: commented these out otherwise the case for multiple MSes would not work
1112 : // Int nsubchans=
1113 : // (chanFreq.shape()(0) - Int(imageStart_p)+1)/Int(imageStep_p)+1;
1114 : // if((nsubchans >0) && (imageNchan_p>nsubchans)) imageNchan_p=nsubchans;
1115 :
1116 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1117 : << "Image spectral coordinate: "<< imageNchan_p
1118 : << " channels, starting at visibility channel "
1119 0 : << imageStart_p << " stepped by " << imageStep_p << LogIO::POST;
1120 0 : freqs.resize(imageNchan_p);
1121 0 : for (Int chan=0;chan<imageNchan_p;chan++) {
1122 0 : freqs(chan)=chanFreq(Int(imageStart_p)+Int(Float(chan)*Float(imageStep_p)));
1123 : }
1124 : // Use this next line when non-linear working
1125 : // mySpectral = new SpectralCoordinate(obsFreqRef, freqs,
1126 : // restFreq);
1127 : // Since we are taking the frequencies as is, the frame must be
1128 : // what is specified in the SPECTRAL_WINDOW table
1129 : // Double finc=(freqs(imageNchan_p-1)-freqs(0))/(imageNchan_p-1);
1130 0 : Double finc = 0;
1131 0 : if(imageNchan_p > 1){
1132 0 : finc=freqs(1)-freqs(0);
1133 : }
1134 0 : else if(imageNchan_p==1) {
1135 0 : finc=freqResolution(IPosition(1,0))*imageStep_p;
1136 : }
1137 :
1138 :
1139 : //in order to outframe to work need to set here original freq frame
1140 : //mySpectral = new SpectralCoordinate(freqFrame_p, freqs(0)-finc/2.0, finc,
1141 0 : mySpectral = new SpectralCoordinate(obsFreqRef, freqs(0)//-finc/2.0
1142 : , finc,
1143 0 : refChan, restFreq);
1144 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1145 : << "Frequency = " // Loglevel INFO
1146 0 : << MFrequency(Quantity(freqs(0), "Hz")).get("GHz").getValue()
1147 : << ", channel increment = "
1148 0 : << MFrequency(Quantity(finc, "Hz")).get("GHz").getValue()
1149 0 : << "GHz" << endl;
1150 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1151 : << "Rest frequency is " // Loglevel INFO
1152 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
1153 0 : << "GHz" << LogIO::POST;
1154 :
1155 0 : }
1156 : // Spectral channels resampled at equal increments in optical velocity
1157 : // Here we compute just the first two channels and use increments for
1158 : // the others
1159 0 : else if (imageMode_p=="VELOCITY" || imageMode_p.contains("RADIO")) {
1160 0 : if(imageNchan_p==0) {
1161 0 : this->unlock();
1162 : os << LogIO::SEVERE << "Must specify number of channels"
1163 0 : << LogIO::EXCEPTION;
1164 0 : return false;
1165 : }
1166 : {
1167 0 : ostringstream oos;
1168 0 : oos << "Image spectral coordinate:"<< imageNchan_p
1169 0 : << " channels, starting at radio velocity " << mImageStart_p
1170 0 : << " stepped by " << mImageStep_p << endl;
1171 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1172 0 : << String(oos); // Loglevel INFO
1173 0 : }
1174 0 : Vector<Double> freqs(2);
1175 0 : freqs=0.0;
1176 0 : if(Double(mImageStep_p.getValue())!=0.0) {
1177 0 : MRadialVelocity mRadVel=mImageStart_p;
1178 0 : for (Int chan=0;chan<2;chan++) {
1179 0 : MDoppler mdoppler(mRadVel.getValue().get(), MDoppler::RADIO);
1180 0 : freqs(chan)=
1181 0 : MFrequency::fromDoppler(mdoppler,
1182 0 : restFreq).getValue().getValue();
1183 0 : Quantity vel=mRadVel.get("m/s");
1184 0 : Quantity inc=mImageStep_p.get("m/s");
1185 0 : vel+=inc;
1186 0 : mRadVel=MRadialVelocity(vel, MRadialVelocity::LSRK);
1187 0 : }
1188 0 : }
1189 : else {
1190 0 : for (Int chan=0;chan<2;++chan) {
1191 0 : freqs(chan)=chanFreq(chan);
1192 : }
1193 : }
1194 :
1195 0 : MFrequency::Types mfreqref=MFrequency::LSRK;
1196 : //Can't convert to frame in mImageStart
1197 0 : if(!MFrequency::getType(mfreqref, (MRadialVelocity::showType(mImageStart_p.getRef().getType()))))
1198 0 : mfreqref=freqFrame_p;
1199 0 : mfreqref=(obsFreqRef==(MFrequency::REST)) ? MFrequency::REST : mfreqref;
1200 0 : mySpectral = new SpectralCoordinate(mfreqref, freqs(0)//-(freqs(1)-freqs(0))/2.0
1201 0 : , freqs(1)-freqs(0), refChan,
1202 0 : restFreq);
1203 :
1204 : {
1205 0 : ostringstream oos;
1206 0 : oos << "Reference Frequency = "
1207 0 : << MFrequency(Quantity(freqs(0), "Hz")).get("GHz")
1208 0 : << ", spectral increment = "
1209 0 : << MFrequency(Quantity(freqs(1)-freqs(0), "Hz")).get("GHz")
1210 0 : << ", frequency frame = "
1211 0 : << MFrequency::showType(mfreqref)
1212 0 : << endl;
1213 0 : oos << "Rest frequency is "
1214 0 : << MFrequency(Quantity(restFreq, "Hz")).get("GHz").getValue()
1215 0 : << " GHz" << endl;
1216 0 : os << LogIO::NORMAL << String(oos) << LogIO::POST; // Loglevel INFO
1217 0 : }
1218 :
1219 0 : }
1220 : // Since optical/relativistic velocity is non-linear in frequency, we have to
1221 : // pass in all the frequencies. For radio velocity we can use
1222 : // a linear axis.
1223 0 : else if (imageMode_p=="OPTICALVELOCITY" || imageMode_p.contains("OPTICAL") || imageMode_p.contains("TRUE")
1224 0 : || imageMode_p.contains("BETA") || imageMode_p.contains("RELATI") ) {
1225 0 : if(imageNchan_p==0) {
1226 0 : this->unlock();
1227 : os << LogIO::SEVERE << "Must specify number of channels"
1228 0 : << LogIO::EXCEPTION;
1229 0 : return false;
1230 : }
1231 : {
1232 0 : ostringstream oos;
1233 0 : oos << "Image spectral coordinate: "<< imageNchan_p
1234 0 : << " channels, starting at optical velocity " << mImageStart_p
1235 0 : << " stepped by " << mImageStep_p << endl;
1236 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1237 0 : << String(oos); // Loglevel INFO
1238 0 : }
1239 0 : Vector<Double> freqs(imageNchan_p);
1240 0 : freqs=0.0;
1241 0 : Double chanVelResolution=0.0;
1242 0 : if(Double(mImageStep_p.getValue())!=0.0) {
1243 0 : MRadialVelocity mRadVel=mImageStart_p;
1244 0 : MDoppler mdoppler;
1245 0 : for (Int chan=0;chan<imageNchan_p;++chan) {
1246 0 : if(imageMode_p.contains("OPTICAL")){
1247 0 : mdoppler=MDoppler(mRadVel.getValue().get(), MDoppler::OPTICAL);
1248 : }
1249 : else{
1250 0 : mdoppler=MDoppler(mRadVel.getValue().get(), MDoppler::BETA);
1251 : }
1252 0 : freqs(chan)=
1253 0 : MFrequency::fromDoppler(mdoppler, restFreq).getValue().getValue();
1254 0 : mRadVel.set(mRadVel.getValue()+mImageStep_p.getValue());
1255 0 : if(imageNchan_p==1)
1256 0 : chanVelResolution=MFrequency::fromDoppler(mdoppler, restFreq).getValue().getValue()-freqs(0);
1257 : }
1258 0 : }
1259 : else {
1260 0 : for (Int chan=0;chan<imageNchan_p;++chan) {
1261 0 : freqs(chan)=chanFreq(chan);
1262 : }
1263 : }
1264 : // Use this next line when non-linear is working
1265 : // when selecting in velocity its specfied freqframe or REST
1266 0 : MFrequency::Types imfreqref=(obsFreqRef==MFrequency::REST) ? MFrequency::REST : freqFrame_p;
1267 :
1268 0 : if(imageNchan_p==1){
1269 0 : if(chanVelResolution==0.0)
1270 0 : chanVelResolution=freqResolution(0);
1271 0 : mySpectral = new SpectralCoordinate(imfreqref,
1272 0 : freqs(0)//-chanVelResolution/2.0,
1273 : ,
1274 : chanVelResolution,
1275 0 : refChan, restFreq);
1276 : }
1277 : else{
1278 0 : mySpectral = new SpectralCoordinate(imfreqref, freqs,
1279 0 : restFreq);
1280 : }
1281 : // mySpectral = new SpectralCoordinate(MFrequency::DEFAULT, freqs(0),
1282 : // freqs(1)-freqs(0), refChan,
1283 : // restFreq);
1284 : {
1285 0 : ostringstream oos;
1286 0 : oos << "Reference Frequency = "
1287 0 : << MFrequency(Quantity(freqs(0), "Hz")).get("GHz")
1288 : << " Ghz, "
1289 0 : <<" frequency frame= "<<MFrequency::showType(imfreqref)<<endl;
1290 : os << (verbose ? LogIO::NORMAL : LogIO::NORMAL3)
1291 0 : << String(oos) << LogIO::POST; // Loglevel INFO
1292 0 : }
1293 0 : }
1294 : else {
1295 0 : this->unlock();
1296 0 : os << LogIO::SEVERE << "Unknown mode " << imageMode_p
1297 0 : << LogIO::EXCEPTION;
1298 0 : return false;
1299 : }
1300 :
1301 :
1302 0 : }
1303 :
1304 :
1305 : //In FTMachine lsrk is used for channel matching with data channel
1306 : //hence we make sure that
1307 : // we convert to lsrk when dealing with the channels
1308 0 : freqFrameValid_p=freqFrameValid_p && (obsFreqRef !=MFrequency::REST);
1309 0 : if(freqFrameValid_p){
1310 0 : mySpectral->setReferenceConversion(MFrequency::LSRK, obsEpoch,
1311 : obsPosition,
1312 0 : phaseCenter_p);
1313 : }
1314 :
1315 : // Polarization
1316 0 : Vector<String> polType=msc.feed().polarizationType()(0);
1317 0 : if (polType(0)!="X" && polType(0)!="Y" &&
1318 0 : polType(0)!="R" && polType(0)!="L") {
1319 : os << LogIO::WARN << "Unknown stokes types in feed table: ["
1320 0 : << polType(0) << ", " << polType(1) << "]" << endl
1321 0 : << "Results open to question!" << LogIO::POST;
1322 : }
1323 :
1324 0 : if (polType(0)=="X" || polType(0)=="Y") {
1325 0 : polRep_p=StokesImageUtil::LINEAR;
1326 : os << LogIO::DEBUG1
1327 0 : << "Preferred polarization representation is linear" << LogIO::POST;
1328 : }
1329 : else {
1330 0 : polRep_p=StokesImageUtil::CIRCULAR;
1331 : os << LogIO::DEBUG1
1332 0 : << "Preferred polarization representation is circular" << LogIO::POST;
1333 : }
1334 :
1335 : // Compare user input with whatever is allowed by the data.
1336 : // If possible, allow.
1337 :
1338 0 : Vector<Int> whichStokes = decideNPolPlanes(true);
1339 0 : if(whichStokes.nelements()==0 || (whichStokes.nelements()==1 && whichStokes[0]==0) )
1340 : {
1341 0 : if(polRep_p==StokesImageUtil::CIRCULAR)
1342 0 : os << LogIO::SEVERE << "Stokes selection of " << stokes_p << " is not valid for Circular feeds." << LogIO::EXCEPTION;
1343 : else
1344 0 : os << LogIO::SEVERE << "Stokes selection of " << stokes_p << " is not valid for Linear feeds." << LogIO::EXCEPTION;
1345 0 : return false;
1346 : }
1347 : /* STOKESDBG */ //cout << "Imager2::imagecoordinate : Image 'whichStokes' : " << whichStokes << endl;
1348 :
1349 0 : StokesCoordinate myStokes(whichStokes);
1350 :
1351 : //Set Observatory info
1352 0 : ObsInfo myobsinfo;
1353 0 : myobsinfo.setTelescope(telescop);
1354 0 : myobsinfo.setPointingCenter(mvPhaseCenter);
1355 0 : myobsinfo.setObsDate(obsEpoch);
1356 0 : myobsinfo.setObserver(msc.observation().observer()(0));
1357 0 : this->setObsInfo(myobsinfo);
1358 :
1359 : //Adding everything to the coordsystem
1360 0 : coordInfo.addCoordinate(myRaDec);
1361 0 : coordInfo.addCoordinate(myStokes);
1362 0 : coordInfo.addCoordinate(*mySpectral);
1363 0 : coordInfo.setObsInfo(myobsinfo);
1364 :
1365 0 : if(mySpectral) delete mySpectral;
1366 :
1367 0 : return true;
1368 0 : }
1369 :
1370 : /* Match the user input stokes_p string to all available options.
1371 : If it matches none, return false - exceptions are thrown.
1372 : If it matches an option, set npol_p : number of desired image planes.
1373 :
1374 : If checkwithMS = true, check with the PolRep and valid combinations of what
1375 : is asked for and what the data allows.
1376 :
1377 : If checkwithMS = false, check only that the inputs are among the valid list
1378 : and give npol=1,2,3 or 4.
1379 :
1380 : There are two levels of checks because when called from defineimage, it does
1381 : not know what PolRep is and cannot check with the data. This used to be in
1382 : two different functions earlier. It is now easier to keep things consistent.
1383 :
1384 : PLEASE NOTE : If any stokes options are to be allowed or disallowed,
1385 : this is the place to do it. StokesImageUtil handles all combinations,
1386 : */
1387 0 : Vector<Int> Imager::decideNPolPlanes(Bool checkwithMS)
1388 : {
1389 0 : Vector<Int> whichStokes(0);
1390 :
1391 : // Check with options for npol=1
1392 0 : if(stokes_p=="I" || stokes_p=="Q" || stokes_p=="U" || stokes_p=="V" ||
1393 0 : stokes_p=="RR" ||stokes_p=="LL" || //stokes_p=="RL" ||stokes_p=="LR" ||
1394 0 : stokes_p=="XX" || stokes_p=="YY" ) //|| stokes_p=="XY" ||stokes_p=="YX" )
1395 : {
1396 0 : npol_p=1;
1397 :
1398 0 : if(checkwithMS)
1399 : {
1400 :
1401 : // Fill in the stokes vector for the output image
1402 0 : whichStokes.resize(npol_p);
1403 : // The first 8 depend on circular vs linear
1404 0 : if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="RR") whichStokes(0)=Stokes::RR;
1405 0 : else if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="LL") whichStokes(0)=Stokes::LL;
1406 : ////else if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="RL") whichStokes(0)=Stokes::RL;
1407 : ////else if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="LR") whichStokes(0)=Stokes::LR;
1408 0 : else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="XX") whichStokes(0)=Stokes::XX;
1409 0 : else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="YY") whichStokes(0)=Stokes::YY;
1410 : ////else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="XY") whichStokes(0)=Stokes::XY;
1411 : ////else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="YX") whichStokes(0)=Stokes::YX;
1412 : // these next 4 don't depend on circular vs linear
1413 0 : else if(stokes_p=="I") whichStokes(0)=Stokes::I;
1414 0 : else if(stokes_p=="Q") whichStokes(0)=Stokes::Q;
1415 0 : else if(stokes_p=="U") whichStokes(0)=Stokes::U;
1416 0 : else if(stokes_p=="V") whichStokes(0)=Stokes::V;
1417 0 : else {whichStokes.resize(0);}
1418 :
1419 : }// end of checkwithMS
1420 :
1421 : }// end of npol=1
1422 : // Check for options with npol=2
1423 0 : else if(stokes_p=="IV" || stokes_p=="IQ" ||
1424 0 : stokes_p=="RRLL" || stokes_p=="XXYY" ||
1425 : ////stokes_p=="RLLR" || stokes_p=="XYYX" ||
1426 0 : stokes_p=="QU" || stokes_p=="UV")
1427 : {
1428 0 : npol_p=2;
1429 :
1430 0 : if(checkwithMS)
1431 : {
1432 : // Check with polrep and fill in the stokes vector for the output image
1433 0 : whichStokes.resize(npol_p);
1434 : // The first 4 depend on circular vs linear
1435 0 : if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="RRLL")
1436 0 : {whichStokes(0)=Stokes::RR; whichStokes(1)=Stokes::LL;}
1437 : ////else if(polRep_p==StokesImageUtil::CIRCULAR && stokes_p=="RLLR")
1438 : //// {whichStokes(0)=Stokes::RL; whichStokes(1)=Stokes::LR;}
1439 0 : else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="XXYY")
1440 0 : {whichStokes(0)=Stokes::XX; whichStokes(1)=Stokes::YY;}
1441 : ////else if(polRep_p==StokesImageUtil::LINEAR && stokes_p=="XYYX")
1442 : //// {whichStokes(0)=Stokes::XY; whichStokes(1)=Stokes::YX;}
1443 : // These 4 don't care about circular vs linear
1444 0 : else if(stokes_p=="IV")
1445 0 : {whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::V;}
1446 0 : else if(stokes_p=="QU")
1447 0 : {whichStokes(0)=Stokes::Q; whichStokes(1)=Stokes::U;}
1448 0 : else if(stokes_p=="IQ")
1449 0 : {whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;}
1450 0 : else if(stokes_p=="UV")
1451 0 : {whichStokes(0)=Stokes::U; whichStokes(1)=Stokes::V;}
1452 : else
1453 0 : {whichStokes.resize(0);}
1454 :
1455 : }// end of checkwithMS
1456 :
1457 : }
1458 : // Check for options with npol=3
1459 0 : else if(stokes_p=="IQU" || stokes_p=="IUV")
1460 : {
1461 0 : npol_p=3;
1462 :
1463 0 : if(checkwithMS)
1464 : {
1465 0 : whichStokes.resize(npol_p);
1466 : // Check with polRep_p
1467 0 : if(stokes_p=="IUV")
1468 0 : { whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::U; whichStokes(2)=Stokes::V; }
1469 0 : else if(stokes_p=="IQU")
1470 0 : { whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q; whichStokes(2)=Stokes::U; }
1471 : else
1472 0 : { whichStokes.resize(0);}
1473 :
1474 : }// end of checkwithMS
1475 :
1476 : }//end of npol=3
1477 : // Check with options for npol=4
1478 0 : else if(stokes_p=="IQUV") // NOT allowing RRLLRLLR and XXYYXYYX !!
1479 : {
1480 0 : npol_p=4;
1481 :
1482 0 : if(checkwithMS)
1483 : {
1484 0 : whichStokes.resize(npol_p);
1485 0 : whichStokes(0)=Stokes::I;
1486 0 : whichStokes(1)=Stokes::Q;
1487 0 : whichStokes(2)=Stokes::U;
1488 0 : whichStokes(3)=Stokes::V;
1489 : }// end of checkwithMS
1490 :
1491 : }// end of npol=4
1492 : else // If none of the options are valid (checked in setimage/defineimage)
1493 : {
1494 0 : whichStokes.resize(1); whichStokes[0]=0; // undefined.
1495 : }
1496 :
1497 : // If no match has been found, or an incompatible match is found, this returns Vector<Int>(0);
1498 0 : return whichStokes;
1499 0 : }// end of decideNPolPlanes
1500 :
1501 :
1502 :
1503 :
1504 0 : String Imager::state()
1505 : {
1506 0 : ostringstream os;
1507 :
1508 : try {
1509 0 : this->lock();
1510 0 : os << "General: " << endl;
1511 0 : os << " MeasurementSet is " << ms_p->tableName() << endl;
1512 0 : if(beamValid_p) {
1513 0 : os << " Beam fit: " << beam_p(0,0).getMajor("arcsec") << " by "
1514 0 : << beam_p(0,0).getMinor("arcsec") << " (arcsec) at pa "
1515 0 : << beam_p(0,0).getPA(Unit("deg")) << " (deg) " << endl;
1516 : }
1517 : else {
1518 0 : os << " Beam fit is not valid" << endl;
1519 : }
1520 :
1521 0 : MSColumns msc(*ms_p);
1522 0 : MDirection mDesiredCenter;
1523 0 : if(setimaged_p) {
1524 : os << "Image definition settings: "
1525 0 : "(use setimage in Function Group <setup> to change)" << endl;
1526 0 : os << " nx=" << nx_p << ", ny=" << ny_p
1527 0 : << ", cellx=" << mcellx_p << ", celly=" << mcelly_p
1528 0 : << ", Stokes axes : " << stokes_p << endl;
1529 0 : Int widthRA=20;
1530 0 : Int widthDec=20;
1531 0 : if(doShift_p) {
1532 0 : os << " doshift is true: Image phase center will be as specified "
1533 0 : << endl;
1534 : }
1535 : else {
1536 0 : os << " doshift is false: Image phase center will be that of field "
1537 0 : << fieldid_p << " :" << endl;
1538 : }
1539 :
1540 0 : if(shiftx_p.get().getValue()!=0.0||shifty_p.get().getValue()!=0.0) {
1541 0 : os << " plus the shift: longitude: " << shiftx_p
1542 0 : << " / cos(latitude) : latitude: " << shifty_p << endl;
1543 : }
1544 :
1545 0 : MVAngle mvRa=phaseCenter_p.getAngle().getValue()(0);
1546 0 : MVAngle mvDec=phaseCenter_p.getAngle().getValue()(1);
1547 0 : os << " ";
1548 0 : os.setf(ios::left, ios::adjustfield);
1549 0 : os.width(widthRA); os << mvRa(0.0).string(MVAngle::TIME,8);
1550 0 : os.width(widthDec); os << mvDec.string(MVAngle::DIG2,8);
1551 0 : os << " " << MDirection::showType(phaseCenter_p.getRefPtr()->getType())
1552 0 : << endl;
1553 :
1554 0 : if(distance_p.get().getValue()!=0.0) {
1555 0 : os << " Refocusing to distance " << distance_p << endl;
1556 : }
1557 :
1558 0 : if(imageMode_p=="MFS") {
1559 0 : os << " Image mode is mfs: Image will be frequency synthesised from spectral windows : ";
1560 0 : for (uInt i=0;i<spectralwindowids_p.nelements();++i) {
1561 0 : os << spectralwindowids_p(i) << " ";
1562 : }
1563 0 : os << endl;
1564 : }
1565 :
1566 : else {
1567 0 : os << " Image mode is " << imageMode_p
1568 0 : << " Image number of spectral channels ="
1569 0 : << imageNchan_p << endl;
1570 : }
1571 :
1572 0 : }
1573 : else {
1574 : os << "Image: parameters not yet set (use setimage "
1575 0 : "in Function Group <setup> )" << endl;
1576 : }
1577 :
1578 : os << "Data selection settings: (use setdata in Function Group <setup> "
1579 0 : "to change)" << endl;
1580 0 : if(dataMode_p=="none") {
1581 0 : if(mssel_p->nrow() == ms_p->nrow()){
1582 0 : os << " All data selected" << endl;
1583 : }
1584 : else{
1585 0 : os << " Number of rows of data selected= " << mssel_p->nrow() << endl;
1586 : }
1587 : }
1588 : else {
1589 0 : os << " Data selection mode is " << dataMode_p << ": have selected "
1590 0 : << dataNchan_p << " channels";
1591 0 : if(dataspectralwindowids_p.nelements()>0) {
1592 0 : os << " spectral windows : ";
1593 0 : for (uInt i=0;i<dataspectralwindowids_p.nelements();++i) {
1594 0 : os << dataspectralwindowids_p(i) << " ";
1595 : }
1596 : }
1597 0 : if(datafieldids_p.nelements()>0) {
1598 0 : os << " Data selected includes fields : ";
1599 0 : for (uInt i=0;i<datafieldids_p.nelements();++i) {
1600 0 : os << datafieldids_p(i) << " ";
1601 : }
1602 : }
1603 0 : os << endl;
1604 : }
1605 : os << "Options settings: (use setoptions in Function Group <setup> "
1606 0 : "to change) " << endl;
1607 0 : os << " Gridding cache has " << cache_p << " complex pixels, in tiles of "
1608 0 : << tile_p << " pixels on a side" << endl;
1609 0 : os << " Gridding convolution function is ";
1610 :
1611 0 : if(gridfunction_p=="SF") {
1612 0 : os << "Spheroidal wave function";
1613 : }
1614 0 : else if(gridfunction_p=="BOX") {
1615 0 : os << "Box car convolution";
1616 : }
1617 0 : else if(gridfunction_p=="PB") {
1618 0 : os << "Using primary beam for convolution";
1619 : }
1620 : else {
1621 0 : os << "Unknown type : " << gridfunction_p;
1622 : }
1623 0 : os << endl;
1624 :
1625 0 : if(doVP_p) {
1626 0 : os << " Primary beam correction is enabled" << endl;
1627 : // Table vpTable( vpTableStr_p ); could fish out info and summarize
1628 : }
1629 : else {
1630 0 : os << " No primary beam correction will be made " << endl;
1631 : }
1632 0 : os << " Image plane padding : " << padding_p << endl;
1633 :
1634 0 : this->unlock();
1635 0 : } catch (AipsError x) {
1636 0 : this->unlock();
1637 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1638 0 : << LogIO::EXCEPTION;
1639 :
1640 0 : }
1641 0 : return String(os);
1642 0 : }
1643 :
1644 :
1645 :
1646 :
1647 : // Apply a primary beam or voltage pattern to an image
1648 0 : Bool Imager::pb(const String& inimage,
1649 : const String& outimage,
1650 : const String& incomps,
1651 : const String& /*outcomps*/,
1652 : const String& operation,
1653 : const MDirection& pointingCenter,
1654 : const Quantity& pa,
1655 : const String& pborvp)
1656 :
1657 : {
1658 0 : if(!valid()) return false;
1659 :
1660 0 : LogIO os(LogOrigin("imager", "pb()", WHERE));
1661 :
1662 0 : PagedImage<Float> * inImage_pointer = 0;
1663 0 : PagedImage<Float> * outImage_pointer = 0;
1664 0 : ComponentList * inComps_pointer = 0;
1665 0 : ComponentList * outComps_pointer = 0;
1666 0 : PBMath * myPBp = 0;
1667 : try {
1668 :
1669 0 : if ( ! doVP_p ) {
1670 0 : this->unlock();
1671 : os << LogIO::SEVERE <<
1672 0 : "Must invoke setvp() first in order to apply the primary beam" << LogIO::EXCEPTION;
1673 0 : return false;
1674 : }
1675 :
1676 0 : if (pborvp == "vp") {
1677 0 : this->unlock();
1678 0 : os << LogIO::SEVERE << "VP application is not yet implemented in DOimager" << LogIO::EXCEPTION;
1679 0 : return false;
1680 : }
1681 :
1682 0 : if (operation == "apply") {
1683 : os << LogIO::DEBUG1
1684 0 : << "function pb will apply " << pborvp << LogIO::POST;
1685 0 : } else if (operation=="correct") {
1686 : os << LogIO::DEBUG1
1687 0 : << "function pb will correct for " << pborvp << LogIO::POST;
1688 : } else {
1689 0 : this->unlock();
1690 : os << LogIO::SEVERE << "Unknown pb operation " << operation
1691 0 : << LogIO::EXCEPTION;
1692 0 : return false;
1693 : }
1694 :
1695 : // Get initial image and/or SkyComponents
1696 :
1697 0 : if (incomps!="") {
1698 0 : if(!Table::isReadable(incomps)) {
1699 0 : this->unlock();
1700 : os << LogIO::SEVERE << "ComponentList " << incomps
1701 0 : << " not readable" << LogIO::EXCEPTION;
1702 0 : return false;
1703 : }
1704 0 : inComps_pointer = new ComponentList(incomps);
1705 : //outComps_pointer = new ComponentList( inComps_pointer->copy() );
1706 : }
1707 0 : if (inimage !="") {
1708 0 : if(!Table::isReadable(inimage)) {
1709 0 : this->unlock();
1710 : os << LogIO::SEVERE << "Image " << inimage << " not readable"
1711 0 : << LogIO::EXCEPTION;
1712 0 : return false;
1713 : }
1714 0 : inImage_pointer = new PagedImage<Float>( inimage );
1715 0 : if (outimage != "") {
1716 0 : outImage_pointer = new PagedImage<Float>( TiledShape(inImage_pointer->shape(), inImage_pointer->niceCursorShape()),
1717 0 : inImage_pointer->coordinates(), outimage);
1718 : }
1719 : }
1720 : // create the PBMath object, needed to apply PB
1721 : // to make high res Fourier weight image
1722 0 : Quantity qFreq;
1723 0 : if (doDefaultVP_p && inImage_pointer==0) {
1724 0 : this->unlock();
1725 : os << LogIO::SEVERE <<
1726 : "There is no default telescope associated with a componentlist"
1727 0 : << LogIO::POST;
1728 : os << LogIO::SEVERE <<
1729 : "Either specify the PB/VP via a vptable or supply an image as well"
1730 0 : << LogIO::EXCEPTION;
1731 0 : return false;
1732 0 : } else if (doDefaultVP_p && inImage_pointer!=0) {
1733 : // look up the telescope in ObsInfo
1734 0 : ObsInfo oi = inImage_pointer->coordinates().obsInfo();
1735 0 : String myTelescope = oi.telescope();
1736 0 : if (myTelescope == "") {
1737 0 : this->unlock();
1738 : os << LogIO::SEVERE << "No telescope imbedded in image"
1739 0 : << LogIO::EXCEPTION;
1740 0 : return false;
1741 : }
1742 : {
1743 0 : Int spectralIndex=inImage_pointer->coordinates().findCoordinate(Coordinate::SPECTRAL);
1744 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1745 : SpectralCoordinate
1746 0 : spectralCoord=inImage_pointer->coordinates().spectralCoordinate(spectralIndex);
1747 0 : Vector<String> units(1); units = "Hz";
1748 0 : spectralCoord.setWorldAxisUnits(units);
1749 0 : Vector<Double> spectralWorld(1);
1750 0 : Vector<Double> spectralPixel(1);
1751 0 : spectralPixel(0) = 0;
1752 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
1753 0 : Double freq = spectralWorld(0);
1754 0 : qFreq = Quantity( freq, "Hz" );
1755 0 : }
1756 0 : String band;
1757 : PBMath::CommonPB whichPB;
1758 0 : String pbName;
1759 : // get freq from coordinates
1760 0 : PBMath::whichCommonPBtoUse (myTelescope, qFreq, band, whichPB, pbName);
1761 0 : if (whichPB == PBMath::UNKNOWN) {
1762 0 : this->unlock();
1763 : os << LogIO::SEVERE << "Unknown telescope for PB type: "
1764 0 : << myTelescope << LogIO::EXCEPTION;
1765 0 : return false;
1766 : }
1767 0 : myPBp = new PBMath(whichPB);
1768 0 : } else {
1769 : // get the PB from the vpTable
1770 0 : Table vpTable( vpTableStr_p );
1771 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
1772 0 : myPBp = new PBMath(recCol(0));
1773 0 : }
1774 0 : AlwaysAssert((myPBp != 0), AipsError);
1775 :
1776 :
1777 : // Do images (if indeed we have any)
1778 0 : if (outImage_pointer!=0) {
1779 0 : Vector<Int> whichStokes;
1780 0 : CoordinateSystem cCoords;
1781 0 : cCoords=StokesImageUtil::CStokesCoord(//inImage_pointer->shape(),
1782 : inImage_pointer->coordinates(),
1783 : whichStokes,
1784 0 : polRep_p);
1785 0 : TempImage<Complex> cIn(inImage_pointer->shape(),
1786 0 : cCoords);
1787 0 : StokesImageUtil::From(cIn, *inImage_pointer);
1788 0 : if (operation=="apply") {
1789 0 : myPBp->applyPB(cIn, cIn, pointingCenter,
1790 : pa, squintType_p, false);
1791 : } else {
1792 0 : myPBp->applyPB(cIn, cIn, pointingCenter,
1793 : pa, squintType_p, true);
1794 : }
1795 0 : StokesImageUtil::To(*outImage_pointer, cIn);
1796 0 : }
1797 : // Do components (if indeed we have any)
1798 0 : if (inComps_pointer!=0) {
1799 0 : if (inImage_pointer==0) {
1800 0 : this->unlock();
1801 : os << LogIO::SEVERE <<
1802 : "No input image was given for the componentList to get the frequency from"
1803 0 : << LogIO::EXCEPTION;
1804 0 : return false;
1805 : }
1806 0 : Int ncomponents = inComps_pointer->nelements();
1807 0 : outComps_pointer = new ComponentList();
1808 0 : for (Int icomp=0;icomp<ncomponents;++icomp) {
1809 0 : SkyComponent component=(inComps_pointer->component(icomp)).copy();
1810 0 : if (operation=="apply") {
1811 0 : myPBp->applyPB(component, component, pointingCenter,
1812 : qFreq, pa, squintType_p, false);
1813 : } else {
1814 0 : myPBp->applyPB(component, component, pointingCenter,
1815 : qFreq, pa, squintType_p, true);
1816 : }
1817 0 : outComps_pointer->add(component);
1818 0 : }
1819 0 : if(outImage_pointer){
1820 0 : ComponentListImage clImage(*outComps_pointer, outImage_pointer->coordinates(), outImage_pointer->shape());
1821 0 : outImage_pointer->copyData(LatticeExpr<Float>((*outImage_pointer)+clImage));
1822 0 : }
1823 : }
1824 0 : if (myPBp) delete myPBp;
1825 0 : if (inImage_pointer) delete inImage_pointer;
1826 0 : if (outImage_pointer) delete outImage_pointer;
1827 0 : if (inComps_pointer) delete inComps_pointer;
1828 0 : if (outComps_pointer) delete outComps_pointer;
1829 0 : return true;
1830 0 : } catch (AipsError x) {
1831 0 : if (myPBp) delete myPBp;
1832 0 : if (inImage_pointer) delete inImage_pointer;
1833 0 : if (outImage_pointer) delete outImage_pointer;
1834 0 : if (inComps_pointer) delete inComps_pointer;
1835 0 : if (outComps_pointer) delete outComps_pointer;
1836 0 : this->unlock();
1837 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1838 0 : << LogIO::EXCEPTION;
1839 0 : return false;
1840 0 : }
1841 : return true;
1842 0 : }
1843 :
1844 :
1845 : // Apply a primary beam or voltage pattern to an image
1846 0 : Bool Imager::pbguts(ImageInterface<Float>& inImage,
1847 : ImageInterface<Float>& outImage,
1848 : const MDirection& pointingDirection,
1849 : const Quantity& pa)
1850 : {
1851 0 : if(!valid()) return false;
1852 :
1853 0 : LogIO os(LogOrigin("imager", "pbguts()", WHERE));
1854 :
1855 : try {
1856 0 : if ( ! doVP_p ) {
1857 : os << LogIO::SEVERE <<
1858 0 : "Must invoke setvp() first in order to apply the primary beam" << LogIO::POST;
1859 0 : return false;
1860 : }
1861 0 : String operation = "apply"; // could have as input in the future!
1862 :
1863 : // create the PBMath object, needed to apply PB
1864 : // to make high res Fourier weight image
1865 0 : Quantity qFreq;
1866 0 : PBMath * myPBp = 0;
1867 :
1868 0 : if (doDefaultVP_p) {
1869 : // look up the telescope in ObsInfo
1870 0 : ObsInfo oi = inImage.coordinates().obsInfo();
1871 0 : String myTelescope = oi.telescope();
1872 0 : if (myTelescope == "") {
1873 0 : os << LogIO::SEVERE << "No telescope embedded in image" << LogIO::POST;
1874 0 : return false;
1875 : }
1876 : {
1877 0 : Int spectralIndex=inImage.coordinates().findCoordinate(Coordinate::SPECTRAL);
1878 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1879 : SpectralCoordinate
1880 0 : spectralCoord=inImage.coordinates().spectralCoordinate(spectralIndex);
1881 0 : Vector<String> units(1); units = "Hz";
1882 0 : spectralCoord.setWorldAxisUnits(units);
1883 0 : Vector<Double> spectralWorld(1);
1884 0 : Vector<Double> spectralPixel(1);
1885 0 : spectralPixel(0) = 0;
1886 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
1887 0 : Double freq = spectralWorld(0);
1888 0 : qFreq = Quantity( freq, "Hz" );
1889 0 : }
1890 0 : String band;
1891 : PBMath::CommonPB whichPB;
1892 0 : String pbName;
1893 : // get freq from coordinates
1894 0 : PBMath::whichCommonPBtoUse (myTelescope, qFreq, band, whichPB, pbName);
1895 0 : if (whichPB != PBMath::UNKNOWN) {
1896 :
1897 0 : myPBp = new PBMath(whichPB);
1898 : }
1899 : else{
1900 0 : MSAntennaColumns ac(ms_p->antenna());
1901 0 : Double dishDiam=ac.dishDiameter()(0);
1902 0 : myPBp= new PBMath(dishDiam, true, qFreq);
1903 :
1904 0 : }
1905 0 : } else {
1906 : // get the PB from the vpTable
1907 0 : Table vpTable( vpTableStr_p );
1908 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
1909 0 : myPBp = new PBMath(recCol(0));
1910 0 : }
1911 0 : AlwaysAssert((myPBp != 0), AipsError);
1912 :
1913 0 : Vector<Int> whichStokes;
1914 0 : CoordinateSystem cCoords;
1915 0 : cCoords=StokesImageUtil::CStokesCoord(//inImage.shape(),
1916 : inImage.coordinates(),
1917 : whichStokes,
1918 0 : polRep_p);
1919 0 : TempImage<Complex> cIn(inImage.shape(),
1920 0 : cCoords);
1921 0 : StokesImageUtil::From(cIn, inImage);
1922 0 : if (operation=="apply") {
1923 0 : myPBp->applyPB(cIn, cIn, pointingDirection,
1924 : pa, squintType_p, false);
1925 : } else {
1926 0 : myPBp->applyPB(cIn, cIn, pointingDirection,
1927 : pa, squintType_p, true);
1928 : }
1929 0 : StokesImageUtil::To(outImage, cIn);
1930 0 : delete myPBp;
1931 0 : return true;
1932 0 : } catch (AipsError x) {
1933 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1934 0 : << LogIO::POST;
1935 0 : return false;
1936 0 : }
1937 :
1938 : return true;
1939 0 : }
1940 :
1941 :
1942 0 : void Imager::printbeam(CleanImageSkyModel *sm_p, LogIO &os, const Bool firstrun)
1943 : {
1944 : //Use predefined beam for restoring or find one by fitting
1945 0 : Bool printBeam = false;
1946 0 : if(!beamValid_p){
1947 0 : ImageBeamSet beam;
1948 0 : beam=sm_p->beam(0);
1949 0 : if(beam.nelements() > 0){
1950 : /*beam_p.setMajorMinor(
1951 : Quantity(abs(beam(0)), "arcsec"),
1952 : Quantity(abs(beam(1)), "arcsec")
1953 : );
1954 : beam_p.setPA(Quantity(beam(2), "deg"));
1955 : */
1956 0 : beam_p=beam;
1957 0 : beamValid_p=true;
1958 0 : printBeam = true;
1959 0 : os << LogIO::NORMAL << "Fitted beam used in restoration: " ; // Loglevel INFO
1960 : }
1961 0 : }
1962 0 : else if(firstrun){
1963 0 : printBeam = true;
1964 0 : os << LogIO::NORMAL << "Beam used in restoration: "; // Loglevel INFO
1965 : }
1966 0 : if(printBeam)
1967 0 : os << LogIO::NORMAL << beam_p(0,0).getMajor("arcsec") << " by " // Loglevel INFO
1968 0 : << beam_p(0,0).getMinor("arcsec") << " (arcsec) at pa "
1969 0 : << beam_p(0,0).getPA(Unit("deg")) << " (deg) " << LogIO::POST;
1970 0 : }
1971 :
1972 0 : Bool Imager::restoreImages(const Vector<String>& restoredNames, Bool modresiduals)
1973 : {
1974 :
1975 0 : LogIO os(LogOrigin("imager", "restoreImages()", WHERE));
1976 : try{
1977 : // It's important that we use the congruent images in both
1978 : // cases. This means that we must use the residual image as
1979 : // passed to the SkyModel and not the one returned. This
1980 : // distinction is important only (currently) for WFCleanImageSkyModel
1981 : // which has a different representation for the image internally.
1982 0 : Vector<String> residualNames(images_p.nelements());
1983 0 : Vector<String> modelNames(images_p.nelements());
1984 0 : Vector<Bool> dofluxscale(images_p.nelements());
1985 0 : dofluxscale=false;
1986 0 : for(uInt k=0; k < modelNames.nelements() ; ++k){
1987 0 : residualNames[k]=residuals_p[k]->name();
1988 0 : modelNames[k]=images_p[k]->name();
1989 0 : dofluxscale[k]=sm_p->doFluxScale(k);
1990 : }
1991 :
1992 :
1993 0 : Double availablemem=Double(HostInfo::memoryFree())*1024.0;
1994 0 : Bool nomemory=false;
1995 0 : Block<CountedPtr< TempImage<Float> > > tempfluximage(modelNames.nelements());
1996 : //The convolution needs in ram 3 complex x-y planes ...lets multiply it by 5 for safety
1997 0 : if((availablemem < Double(nx_p*ny_p)*15.0*8.0 ) && (ft_p->name() != "MosaicFT") && !(doWideBand_p && ntaylor_p>1)){
1998 : // very large for convolution ...destroy Skyequations to release memory
1999 : // need to fix the convolution to be leaner
2000 0 : for (Int thismodel=0;thismodel<Int(restoredNames.nelements());
2001 : ++thismodel) {
2002 0 : tempfluximage[thismodel]=new TempImage<Float>(sm_p->fluxScale(thismodel).shape(), sm_p->fluxScale(thismodel).coordinates(), 10.0);
2003 0 : (tempfluximage[thismodel])->copyData(sm_p->fluxScale(thismodel));
2004 : }
2005 0 : destroySkyEquation();
2006 0 : nomemory=true;
2007 :
2008 : }
2009 :
2010 :
2011 : // If msmfs, calculate Coeff Residuals
2012 0 : if(doWideBand_p && ntaylor_p>1)
2013 : {
2014 0 : if( modresiduals ) // When called from pclean, this is set to false, via the iClean call to restoreImages.
2015 : {
2016 0 : sm_p->calculateCoeffResiduals();
2017 : // Re-fill them into the output residual images.
2018 0 : for (uInt k=0 ; k < residuals_p.nelements(); ++k){
2019 0 : (residuals_p[k])->copyData(sm_p->getResidual(k));
2020 : }
2021 : }
2022 : }
2023 :
2024 :
2025 0 : Bool dorestore=false;
2026 0 : if( beam_p.nelements() >0 )
2027 0 : dorestore=true;
2028 :
2029 0 : if(restoredNames.nelements()>0) {
2030 0 : for (Int thismodel=0;thismodel<Int(restoredNames.nelements());
2031 : ++thismodel) {
2032 0 : if(restoredNames(thismodel)!="") {
2033 0 : PagedImage<Float> modelIm(modelNames[thismodel]);
2034 0 : PagedImage<Float> residIm(residualNames[thismodel]);
2035 0 : TempImage<Float> restored;
2036 0 : if(dorestore){
2037 0 : restored=TempImage<Float>(modelIm.shape(),
2038 0 : modelIm.coordinates());
2039 0 : restored.copyData(modelIm);
2040 :
2041 0 : StokesImageUtil::Convolve(restored, beam_p);
2042 : }
2043 : // We can work only if the residual image was defined.
2044 0 : if(residIm.name() != "") {
2045 0 : if(dorestore){
2046 0 : LatticeExpr<Float> le(restored+(residIm));
2047 0 : restored.copyData(le);
2048 0 : }
2049 0 : if(freqFrameValid_p){
2050 0 : CoordinateSystem cs=residIm.coordinates();
2051 0 : String errorMsg;
2052 0 : if (cs.setSpectralConversion (errorMsg, MFrequency::showType(freqFrame_p))) {
2053 0 : residIm.setCoordinateInfo(cs);
2054 0 : if(dorestore)
2055 0 : restored.setCoordinateInfo(cs);
2056 : }
2057 0 : }
2058 :
2059 : //should be able to do that only on testing dofluxscale
2060 : // ftmachines or sm_p should tell us that
2061 :
2062 : // this should go away..
2063 : //special casing like this gets hard to maintain
2064 : // need to redo how interactive is done so that it is outside
2065 : //of this code
2066 :
2067 :
2068 : //
2069 : // Using minPB_p^2 below to make it consistent with the normalization in SkyEquation for ftmachine mosaic as coverimage in that case in square of what it should be
2070 : //
2071 0 : Float cutoffval=minPB_p;
2072 0 : if(ft_p->name()=="MosaicFT")
2073 0 : cutoffval=minPB_p*minPB_p;
2074 :
2075 0 : ImageInterface<Float> * diviseur= nomemory ? &(*tempfluximage[thismodel]) : &(sm_p->fluxScale(thismodel));
2076 :
2077 0 : if (dofluxscale(thismodel)) {
2078 0 : TempImage<Float> cover(modelIm.shape(),modelIm.coordinates());
2079 0 : if(ft_p->name()=="MosaicFT")
2080 0 : se_p->getCoverageImage(thismodel, cover);
2081 : else
2082 0 : cover.copyData(*diviseur);
2083 :
2084 0 : if(scaleType_p=="NONE"){
2085 0 : if(dorestore){
2086 0 : LatticeExpr<Float> le(iif(((cover > cutoffval) && ((*diviseur) >0.0)) ,
2087 0 : (restored/(*diviseur)), 0.0));
2088 0 : restored.copyData(le);
2089 0 : }
2090 0 : LatticeExpr<Float> le1(iif(((cover > cutoffval) && ((*diviseur) >0.0)),
2091 0 : (residIm/(*diviseur)), 0.0));
2092 0 : residIm.copyData(le1);
2093 :
2094 0 : }
2095 :
2096 : //Setting the bit-mask for mosaic image
2097 0 : LatticeExpr<Bool> lemask(iif((cover > cutoffval) && ((*diviseur) >0.0),
2098 0 : true, false));
2099 0 : if(dorestore){
2100 0 : ImageRegion outreg=restored.makeMask("mask0", false, true);
2101 0 : LCRegion& outmask=outreg.asMask();
2102 0 : outmask.copyData(lemask);
2103 0 : restored.defineRegion("mask0", outreg, RegionHandler::Masks, true);
2104 0 : restored.setDefaultMask("mask0");
2105 :
2106 0 : }
2107 :
2108 0 : }
2109 0 : if(dorestore){
2110 0 : PagedImage<Float> diskrestore(restored.shape(), restored.coordinates(),
2111 0 : restoredNames(thismodel));
2112 0 : diskrestore.copyData(restored);
2113 0 : ImageInfo ii = modelIm.imageInfo();
2114 0 : ii.setBeams(beam_p);
2115 0 : diskrestore.setImageInfo(ii);
2116 0 : diskrestore.setUnits(Unit("Jy/beam"));
2117 0 : copyMask(diskrestore, restored, "mask0");
2118 :
2119 0 : }
2120 :
2121 : }
2122 : else {
2123 : os << LogIO::SEVERE << "No residual image for model "
2124 0 : << thismodel << ", cannot restore image" << LogIO::POST;
2125 : }
2126 :
2127 0 : if(!(residuals_p[thismodel].null()) && residuals_p[thismodel]->ok()){
2128 0 : residuals_p[thismodel]->table().relinquishAutoLocks(true);
2129 0 : residuals_p[thismodel]->table().unlock();
2130 : //detaching residual so that other processes can use it
2131 0 : residuals_p[thismodel]=0;
2132 : }
2133 0 : }
2134 :
2135 : }// end of for 'thismodel'
2136 :
2137 : // If msmfs, calculate alpha, beta too
2138 0 : if(doWideBand_p && ntaylor_p>1)
2139 : {
2140 0 : sm_p->calculateAlphaBeta(restoredNames, residualNames);
2141 : }
2142 :
2143 : }
2144 0 : }
2145 0 : catch (exception &x) {
2146 0 : os << LogIO::SEVERE << "Exception: " << x.what() << LogIO::POST;
2147 0 : return false;
2148 0 : }
2149 0 : return true;
2150 0 : }
2151 :
2152 0 : Bool Imager::copyMask(ImageInterface<Float>& out, const ImageInterface<Float>& in, String maskname, Bool setdef){
2153 : try{
2154 0 : if(in.hasRegion(maskname) && !out.hasRegion(maskname)){
2155 0 : ImageRegion outreg=out.makeMask(maskname, false, true);
2156 0 : LCRegion& outmask=outreg.asMask();
2157 0 : outmask.copyData(in.getRegion(maskname).asLCRegion());
2158 0 : LatticeExpr<Float> myexpr(iif(outmask, out, 0.0) );
2159 0 : out.copyData(myexpr);
2160 0 : out.defineRegion(maskname, outreg, RegionHandler::Masks, true);
2161 0 : if(setdef)
2162 0 : out.setDefaultMask(maskname);
2163 0 : }
2164 : else
2165 0 : return false;
2166 :
2167 : }
2168 0 : catch(exception &x){
2169 0 : throw(AipsError(x.what()));
2170 : return false;
2171 0 : }
2172 :
2173 :
2174 0 : return true;
2175 : }
2176 :
2177 0 : Bool Imager::writeFluxScales(const Vector<String>& fluxScaleNames)
2178 : {
2179 0 : LogIO os(LogOrigin("imager", "writeFluxScales()", WHERE));
2180 0 : Bool answer = false;
2181 : ImageInterface<Float> *cover;
2182 0 : if(fluxScaleNames.nelements()>0) {
2183 0 : for (Int thismodel=0;thismodel<Int(fluxScaleNames.nelements());++thismodel) {
2184 0 : if(fluxScaleNames(thismodel)!="") {
2185 0 : PagedImage<Float> fluxScale(images_p[thismodel]->shape(),
2186 0 : images_p[thismodel]->coordinates(),
2187 0 : fluxScaleNames(thismodel));
2188 0 : PagedImage<Float> coverimage(images_p[thismodel]->shape(),
2189 0 : images_p[thismodel]->coordinates(),
2190 0 : fluxScaleNames(thismodel)+".pbcoverage");
2191 0 : coverimage.table().markForDelete();
2192 0 : if(freqFrameValid_p){
2193 0 : CoordinateSystem cs=fluxScale.coordinates();
2194 0 : String errorMsg;
2195 0 : if (cs.setSpectralConversion (errorMsg,MFrequency::showType(freqFrame_p))) {
2196 0 : fluxScale.setCoordinateInfo(cs);
2197 0 : coverimage.setCoordinateInfo(cs);
2198 : }
2199 0 : }
2200 0 : if (sm_p->doFluxScale(thismodel)) {
2201 0 : cover=&(sm_p->fluxScale(thismodel));
2202 0 : answer = true;
2203 0 : fluxScale.copyData(sm_p->fluxScale(thismodel));
2204 0 : Float cutoffval=minPB_p;
2205 0 : if(ft_p->name()=="MosaicFT"){
2206 0 : cutoffval=minPB_p*minPB_p;
2207 0 : se_p->getCoverageImage(thismodel, coverimage);
2208 0 : cover=&(coverimage);
2209 : //Do the sqrt
2210 0 : coverimage.copyData(( LatticeExpr<Float> )(iif(coverimage > 0.0, sqrt(coverimage), 0.0)));
2211 0 : coverimage.table().unmarkForDelete();
2212 0 : LatticeExpr<Bool> lemask(iif((*cover) < sqrt(cutoffval),
2213 0 : false, true));
2214 0 : ImageRegion outreg=coverimage.makeMask("mask0", false, true);
2215 0 : LCRegion& outmask=outreg.asMask();
2216 0 : outmask.copyData(lemask);
2217 0 : coverimage.defineRegion("mask0", outreg,RegionHandler::Masks, true);
2218 0 : coverimage.setDefaultMask("mask0");
2219 0 : }
2220 0 : LatticeExpr<Bool> lemask(iif(((*cover) > minPB_p) && (fluxScale > 0.0),
2221 0 : true, false));
2222 0 : ImageRegion outreg=fluxScale.makeMask("mask0", false, true);
2223 0 : LCRegion& outmask=outreg.asMask();
2224 0 : outmask.copyData(lemask);
2225 0 : fluxScale.defineRegion("mask0", outreg,RegionHandler::Masks, true);
2226 0 : fluxScale.setDefaultMask("mask0");
2227 :
2228 :
2229 0 : } else {
2230 0 : answer = false;
2231 : os << LogIO::NORMAL // Loglevel INFO
2232 : << "No flux scale available (or required) for model " << thismodel
2233 0 : << LogIO::POST;
2234 : os << LogIO::NORMAL // Loglevel INFO
2235 0 : << "(This is only pertinent to mosaiced images)" << LogIO::POST;
2236 : os << LogIO::NORMAL // Loglevel INFO
2237 0 : << "Writing out image of constant 1.0" << LogIO::POST;
2238 0 : fluxScale.set(1.0);
2239 : }
2240 0 : }
2241 : }
2242 : }
2243 0 : return answer;
2244 0 : }
2245 : // Supports the "[] or -1 => everything" convention using the rule:
2246 : // If v is empty or only has 1 element, and it is < 0,
2247 : // replace v with 0, 1, ..., nelem - 1.
2248 : // Returns whether or not it modified v.
2249 : // If so, v is modified in place.
2250 14 : Bool Imager::expand_blank_sel(Vector<Int>& v, const uInt nelem)
2251 : {
2252 14 : if((v.nelements() == 1 && v[0] < 0) || v.nelements() == 0){
2253 7 : v.resize(nelem);
2254 7 : indgen(v);
2255 7 : return true;
2256 : }
2257 7 : return false;
2258 : }
2259 :
2260 :
2261 : // Correct the data using a plain VisEquation.
2262 : // This just moves data from observed to corrected.
2263 : // Eventually we should pass in a calibrater
2264 : // object to do the work.
2265 0 : Bool Imager::correct(const Bool /*doparallactic*/, const Quantity& /*t*/)
2266 : {
2267 0 : if(!valid()) return false;
2268 :
2269 0 : LogIO os(LogOrigin("imager", "correct()", WHERE));
2270 :
2271 : // (gmoellen 06Nov20)
2272 : // VisEquation/VisJones have been extensively revised, so
2273 : // so we are disabling this method. Will probably soon
2274 : /// delete it entirely, since this is a calibrater responsibility....
2275 :
2276 0 : this->lock();
2277 : try {
2278 :
2279 : // Warn users we are disabled.
2280 0 : throw(AipsError("Imager::correct() has been disabled (gmoellen 06Nov20)."));
2281 :
2282 : /* Commenting out old VisEquation/VisJones usage
2283 : os << "Correcting data: CORRECTED_DATA column will be replaced"
2284 : << LogIO::POST;
2285 :
2286 :
2287 : if(doparallactic) {
2288 : os<<"Correcting parallactic angle variation"<<LogIO::POST;
2289 : VisEquation ve(*vs_p);
2290 : Float ts=t.get("s").getValue();
2291 : Float dt=ts/10.0;
2292 : PJones pj(*vs_p, ts, dt);
2293 : ve.setVisJones(pj);
2294 : ve.correct();
2295 : }
2296 : else {
2297 : VisEquation ve(*vs_p);
2298 : ve.correct();
2299 : }
2300 : this->unlock();
2301 : return true;
2302 : */
2303 :
2304 :
2305 0 : } catch (AipsError x) {
2306 0 : this->unlock();
2307 0 : os << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
2308 0 : return false;
2309 0 : }
2310 : this->unlock();
2311 :
2312 : return true;
2313 0 : }
2314 :
2315 0 : uInt Imager::count_visibilities(ROVisibilityIterator *rvi,
2316 : const Bool unflagged_only,
2317 : const Bool must_have_imwt)
2318 : {
2319 0 : if(!valid()) return 0;
2320 0 : LogIO os(LogOrigin("imager", "count_visibilities()", WHERE));
2321 :
2322 0 : this->lock();
2323 0 : ROVisIter& vi(*rvi);
2324 0 : VisBuffer vb(vi);
2325 :
2326 0 : uInt nVis = 0;
2327 0 : for(vi.originChunks(); vi.moreChunks(); vi.nextChunk()){
2328 0 : for(vi.origin(); vi.more(); ++vi){
2329 0 : uInt nRow = vb.nRow();
2330 0 : uInt nChan = vb.nChannel();
2331 0 : for(uInt row = 0; row < nRow; ++row)
2332 0 : for(uInt chn = 0; chn < nChan; ++chn)
2333 0 : if((!unflagged_only || !vb.flag()(chn,row)) &&
2334 0 : (!must_have_imwt || vb.imagingWeight()(chn,row) > 0.0))
2335 0 : ++nVis;
2336 : }
2337 : }
2338 0 : this->unlock();
2339 0 : return nVis;
2340 0 : }
2341 : // Create the FTMachine as late as possible
2342 37 : Bool Imager::createFTMachine()
2343 : {
2344 :
2345 :
2346 37 : if(ft_p) {delete ft_p; ft_p=0;}
2347 37 : if(gvp_p) {delete gvp_p; gvp_p=0;}
2348 37 : if(cft_p) {delete cft_p; cft_p=0;}
2349 :
2350 : //For ftmachines that can use double precision gridders
2351 37 : Bool useDoublePrecGrid=false;
2352 : //few channels use Double precision
2353 : //till we find a better algorithm to determine when to use Double prec gridding
2354 37 : if((imageNchan_p < 5) && !(singlePrec_p))
2355 37 : useDoublePrecGrid=true;
2356 :
2357 74 : LogIO os(LogOrigin("imager", "createFTMachine()", WHERE));
2358 :
2359 : Float padding;
2360 37 : padding=1.0;
2361 37 : if(doMultiFields_p||(facets_p>1)) {
2362 0 : padding = padding_p;
2363 : os << LogIO::NORMAL // Loglevel INFO
2364 : << "Multiple fields or facets: transforms will be padded by a factor "
2365 0 : << padding << LogIO::POST;
2366 : }
2367 :
2368 37 : if(ftmachine_p=="sd") {
2369 : os << LogIO::NORMAL // Loglevel INFO
2370 0 : << "Performing single dish gridding..." << LogIO::POST;
2371 : os << LogIO::NORMAL1 // gridfunction_p is too cryptic for most users.
2372 0 : << "with convolution function " << gridfunction_p << LogIO::POST;
2373 :
2374 : // Now make the Single Dish Gridding
2375 : os << LogIO::NORMAL // Loglevel INFO
2376 0 : << "Gridding will use specified common tangent point:" << LogIO::POST;
2377 0 : os << LogIO::NORMAL << tangentPoint() << LogIO::POST; // Loglevel INFO
2378 0 : if(gridfunction_p=="pb") {
2379 0 : if(!gvp_p) {
2380 0 : MSColumns msc(*ms_p);
2381 0 : if (doDefaultVP_p) {
2382 : os << LogIO::NORMAL // Loglevel INFO
2383 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2384 0 : gvp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p,
2385 0 : skyPosThreshold_p);
2386 : } else {
2387 : os << LogIO::NORMAL // Loglevel INFO
2388 0 : << "Using VP as defined in " << vpTableStr_p << LogIO::POST;
2389 0 : Table vpTable( vpTableStr_p );
2390 0 : gvp_p=new VPSkyJones(msc, vpTable, parAngleInc_p, squintType_p,
2391 0 : skyPosThreshold_p);
2392 0 : }
2393 0 : }
2394 0 : ft_p = new SDGrid(mLocation_p, *gvp_p, cache_p/2, tile_p, gridfunction_p,
2395 0 : sdConvSupport_p, minWeight_p, clipminmax_p);
2396 : }
2397 0 : else if (gridfunction_p=="gauss" || gridfunction_p=="gjinc") {
2398 0 : if (mcellx_p != mcelly_p &&
2399 0 : ((!qtruncate_p.getUnit().empty()||qtruncate_p.getUnit()=="pixel")
2400 0 : || (!qgwidth_p.getUnit().empty()||qgwidth_p.getUnit()=="pixel")
2401 0 : || (!qjwidth_p.getUnit().empty()||qjwidth_p.getUnit()=="pixel"))) {
2402 : os << LogIO::WARN
2403 0 : << "The " << gridfunction_p << " gridding doesn't support non-square grid." << endl
2404 0 : << "Result may be wrong." << LogIO::POST;
2405 : }
2406 : Float truncate, gwidth, jwidth;
2407 0 : if (qtruncate_p.getUnit().empty() || qtruncate_p.getUnit()=="pixel")
2408 0 : truncate = qtruncate_p.getValue();
2409 : else
2410 0 : truncate = qtruncate_p.getValue("rad")/mcelly_p.getValue("rad");
2411 0 : if (qgwidth_p.getUnit().empty() || qgwidth_p.getUnit()=="pixel")
2412 0 : gwidth = qgwidth_p.getValue();
2413 : else
2414 0 : gwidth = qgwidth_p.getValue("rad")/mcelly_p.getValue("rad");
2415 0 : if (qjwidth_p.getUnit().empty() || qjwidth_p.getUnit()=="pixel")
2416 0 : jwidth = qjwidth_p.getValue();
2417 : else
2418 0 : jwidth = qjwidth_p.getValue("rad")/mcelly_p.getValue("rad");
2419 0 : ft_p = new SDGrid(mLocation_p, cache_p/2, tile_p, gridfunction_p,
2420 0 : truncate, gwidth, jwidth, minWeight_p, clipminmax_p);
2421 : }
2422 : else {
2423 0 : ft_p = new SDGrid(mLocation_p, cache_p/2, tile_p, gridfunction_p,
2424 0 : sdConvSupport_p, minWeight_p, clipminmax_p);
2425 : }
2426 0 : ft_p->setPointingDirColumn(pointingDirCol_p);
2427 0 : auto * sdgrid = dynamic_cast<SDGrid *>(ft_p);
2428 0 : if (sdgrid) {
2429 0 : sdgrid->setEnableCache(enablecache_p);
2430 0 : sdgrid->setConvertFirst(convertfirst_p);
2431 : }
2432 :
2433 0 : ROVisIter& vi(*rvi_p);
2434 : // Get bigger chunks o'data: this should be tuned some time
2435 : // since it may be wrong for e.g. spectral line
2436 0 : vi.setRowBlocking(100);
2437 :
2438 0 : AlwaysAssert(ft_p, AipsError);
2439 :
2440 0 : cft_p = new SimpleComponentGridMachine();
2441 0 : AlwaysAssert(cft_p, AipsError);
2442 : }
2443 37 : else if(ftmachine_p=="mosaic") {
2444 0 : os << LogIO::NORMAL << "Performing mosaic gridding" << LogIO::POST; // Loglevel PROGRESS
2445 :
2446 0 : setMosaicFTMachine(useDoublePrecGrid);
2447 :
2448 : // VisIter& vi(vs_p->iter());
2449 : // vi.setRowBlocking(100);
2450 :
2451 0 : AlwaysAssert(ft_p, AipsError);
2452 :
2453 0 : cft_p = new SimpleComponentFTMachine();
2454 0 : AlwaysAssert(cft_p, AipsError);
2455 : }
2456 : //
2457 : // Make WProject FT machine (for non co-planar imaging)
2458 : //
2459 37 : else if (ftmachine_p == "wproject"){
2460 : os << LogIO::NORMAL << "Performing w-plane projection" // Loglevel PROGRESS
2461 0 : << LogIO::POST;
2462 :
2463 0 : Double maxW=-1.0;
2464 0 : Double minW=-1.0;
2465 0 : Double rmsW=-1.0;
2466 0 : if(wprojPlanes_p < 1)
2467 0 : WProjectFT::wStat(*rvi_p, minW, maxW, rmsW);
2468 0 : if(facets_p > 1){
2469 0 : ft_p = new WProjectFT(wprojPlanes_p, phaseCenter_p, mLocation_p,
2470 0 : cache_p/2, tile_p, false, padding_p, useDoublePrecGrid, minW, maxW, rmsW);
2471 : }
2472 : else{
2473 0 : ft_p = new WProjectFT(wprojPlanes_p, mLocation_p,
2474 0 : cache_p/2, tile_p, false, padding_p, useDoublePrecGrid, minW, maxW, rmsW);
2475 : }
2476 0 : AlwaysAssert(ft_p, AipsError);
2477 0 : cft_p = new SimpleComponentFTMachine();
2478 0 : AlwaysAssert(cft_p, AipsError);
2479 : }
2480 : //
2481 : // Make PBWProject FT machine (for non co-planar imaging with
2482 : // antenna based PB corrections)
2483 : //
2484 37 : else if (ftmachine_p == "pbwproject"){
2485 0 : if (wprojPlanes_p<=1)
2486 : {
2487 : os << LogIO::NORMAL
2488 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
2489 0 : << LogIO::POST;
2490 0 : os << LogIO::NORMAL << "Performing pb-projection" << LogIO::POST; // Loglevel PROGRESS
2491 : }
2492 0 : if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
2493 : {
2494 : os << LogIO::WARN
2495 : << "No. of w-planes set too low for W projection - recommend at least 128"
2496 0 : << LogIO::POST;
2497 : os << LogIO::NORMAL << "Performing pb + w-plane projection" // Loglevel PROGRESS
2498 0 : << LogIO::POST;
2499 : }
2500 :
2501 :
2502 0 : if(!gvp_p)
2503 : {
2504 : os << LogIO::NORMAL // Loglevel INFO
2505 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2506 0 : MSColumns msc(*ms_p);
2507 0 : gvp_p = new VPSkyJones(msc, true, parAngleInc_p, squintType_p,
2508 0 : skyPosThreshold_p);
2509 0 : }
2510 :
2511 0 : ft_p = new nPBWProjectFT(//*ms_p,
2512 0 : wprojPlanes_p, cache_p/2,
2513 0 : cfCacheDirName_p, doPointing, doPBCorr,
2514 0 : tile_p, computePAStep_p, pbLimit_p, true);
2515 : //
2516 : // Explicit type casting of ft_p does not look good. It does not
2517 : // pick up the setPAIncrement() method of PBWProjectFT without
2518 : // this
2519 : //
2520 0 : ((nPBWProjectFT *)ft_p)->setPAIncrement(parAngleInc_p);
2521 0 : if (doPointing)
2522 : {
2523 : try
2524 : {
2525 : // epJ = new EPJones(*vs_p, *ms_p);
2526 0 : VisSet elVS(*rvi_p);
2527 0 : epJ = new EPJones(elVS);
2528 0 : RecordDesc applyRecDesc;
2529 0 : applyRecDesc.addField("table", TpString);
2530 0 : applyRecDesc.addField("interp",TpString);
2531 0 : Record applyRec(applyRecDesc);
2532 0 : applyRec.define("table",epJTableName_p);
2533 0 : applyRec.define("interp", "nearest");
2534 0 : epJ->setApply(applyRec);
2535 0 : ((nPBWProjectFT *)ft_p)->setEPJones(epJ);
2536 0 : }
2537 0 : catch(AipsError& x)
2538 : {
2539 : //
2540 : // Add some more useful info. to the message and translate
2541 : // the generic AipsError exception object to a more specific
2542 : // SynthesisError object.
2543 : //
2544 0 : String mesg = x.getMesg();
2545 0 : mesg += ". Error in loading pointing offset table.";
2546 0 : SynthesisError err(mesg);
2547 0 : throw(err);
2548 0 : }
2549 : }
2550 :
2551 0 : AlwaysAssert(ft_p, AipsError);
2552 0 : cft_p = new SimpleComponentFTMachine();
2553 0 : AlwaysAssert(cft_p, AipsError);
2554 : }
2555 37 : else if (ftmachine_p == "pbmosaic"){
2556 :
2557 0 : if (wprojPlanes_p<=1)
2558 : {
2559 : os << LogIO::NORMAL
2560 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
2561 0 : << LogIO::POST;
2562 0 : os << LogIO::NORMAL << "Performing pb-mosaic" << LogIO::POST; // Loglevel PROGRESS
2563 : }
2564 0 : if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
2565 : {
2566 : os << LogIO::WARN
2567 : << "No. of w-planes set too low for W projection - recommend at least 128"
2568 0 : << LogIO::POST;
2569 0 : os << LogIO::NORMAL << "Performing pb + w-plane projection" << LogIO::POST; // Loglevel PROGRESS
2570 : }
2571 :
2572 0 : if(!gvp_p)
2573 : {
2574 : os << LogIO::NORMAL // Loglevel INFO
2575 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2576 0 : MSColumns msc(*ms_p);
2577 0 : gvp_p = new VPSkyJones(msc, true, parAngleInc_p, squintType_p,
2578 0 : skyPosThreshold_p);
2579 0 : }
2580 0 : ft_p = new PBMosaicFT(*ms_p, wprojPlanes_p, cache_p/2,
2581 0 : cfCacheDirName_p, /*true */doPointing, doPBCorr,
2582 0 : tile_p, computePAStep_p, pbLimit_p, true);
2583 0 : ((PBMosaicFT *)ft_p)->setObservatoryLocation(mLocation_p);
2584 : //
2585 : // Explicit type casting of ft_p does not look good. It does not
2586 : // pick up the setPAIncrement() method of PBWProjectFT without
2587 : // this
2588 : //
2589 0 : os << LogIO::NORMAL << "Setting PA increment to " << parAngleInc_p.getValue("deg") << " deg" << endl;
2590 0 : ((nPBWProjectFT *)ft_p)->setPAIncrement(parAngleInc_p);
2591 :
2592 0 : if (doPointing)
2593 : {
2594 : try
2595 : {
2596 : // Warn users we are have temporarily disabled pointing cal
2597 : // throw(AipsError("Pointing calibration temporarily disabled (gmoellen 06Nov20)."));
2598 : // TBD: Bring this up-to-date with new VisCal mechanisms
2599 0 : VisSet elVS(*rvi_p);
2600 0 : epJ = new EPJones(elVS, *ms_p);
2601 0 : RecordDesc applyRecDesc;
2602 0 : applyRecDesc.addField("table", TpString);
2603 0 : applyRecDesc.addField("interp",TpString);
2604 0 : Record applyRec(applyRecDesc);
2605 0 : applyRec.define("table",epJTableName_p);
2606 0 : applyRec.define("interp", "nearest");
2607 0 : epJ->setApply(applyRec);
2608 0 : ((nPBWProjectFT *)ft_p)->setEPJones(epJ);
2609 0 : }
2610 0 : catch(AipsError& x)
2611 : {
2612 : //
2613 : // Add some more useful info. to the message and translate
2614 : // the generic AipsError exception object to a more specific
2615 : // SynthesisError object.
2616 : //
2617 0 : String mesg = x.getMesg();
2618 0 : mesg += ". Error in loading pointing offset table.";
2619 0 : SynthesisError err(mesg);
2620 0 : throw(err);
2621 0 : }
2622 : }
2623 :
2624 0 : AlwaysAssert(ft_p, AipsError);
2625 0 : cft_p = new SimpleComponentFTMachine();
2626 0 : AlwaysAssert(cft_p, AipsError);
2627 :
2628 : }
2629 37 : else if(ftmachine_p=="both") {
2630 :
2631 : os << LogIO::NORMAL // Loglevel INFO
2632 : << "Performing single dish gridding with convolution function "
2633 0 : << gridfunction_p << LogIO::POST;
2634 : os << LogIO::NORMAL // Loglevel INFO
2635 : << "and interferometric gridding with the prolate spheroidal convolution function"
2636 0 : << LogIO::POST;
2637 :
2638 : // Now make the Single Dish Gridding
2639 : os << LogIO::NORMAL // Loglevel INFO
2640 0 : << "Gridding will use specified common tangent point:" << LogIO::POST;
2641 0 : os << LogIO::NORMAL << tangentPoint() << LogIO::POST; // Loglevel INFO
2642 0 : if(!gvp_p) {
2643 : os << LogIO::NORMAL // Loglevel INFO
2644 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2645 0 : MSColumns msc(*ms_p);
2646 0 : gvp_p = new VPSkyJones(msc, true, parAngleInc_p, squintType_p,
2647 0 : skyPosThreshold_p);
2648 0 : }
2649 0 : if(sdScale_p != 1.0)
2650 : os << LogIO::NORMAL // Loglevel INFO
2651 0 : << "Multiplying single dish data by factor " << sdScale_p << LogIO::POST;
2652 0 : if(sdWeight_p != 1.0)
2653 : os << LogIO::NORMAL // Loglevel INFO
2654 0 : << "Multiplying single dish weights by factor " << sdWeight_p
2655 0 : << LogIO::POST;
2656 0 : ft_p = new GridBoth(*gvp_p, cache_p/2, tile_p,
2657 0 : mLocation_p, phaseCenter_p,
2658 0 : gridfunction_p, "SF", padding,
2659 0 : sdScale_p, sdWeight_p);
2660 :
2661 : //VisIter& vi(vs_p->iter());
2662 : // Get bigger chunks o'data: this should be tuned some time
2663 : // since it may be wrong for e.g. spectral line
2664 0 : rvi_p->setRowBlocking(100);
2665 :
2666 0 : AlwaysAssert(ft_p, AipsError);
2667 :
2668 0 : cft_p = new SimpleComponentFTMachine();
2669 0 : AlwaysAssert(cft_p, AipsError);
2670 :
2671 : }
2672 37 : else if(ftmachine_p=="nift") {
2673 : os << LogIO::NORMAL // Loglevel INFO
2674 0 : << "Using FTMachine " << ftmachine_p << LogIO::POST
2675 : << "Performing interferometric gridding..."
2676 0 : << LogIO::POST;
2677 : os << LogIO::NORMAL1 // gridfunction_p is too cryptic for most users.
2678 0 : << "...with convolution function " << gridfunction_p << LogIO::POST;
2679 : /*
2680 : // Make the re-gridder components. Here, make the basic
2681 : // re-sampler.
2682 : CountedPtr<VisibilityResamplerBase> visResamplerCtor = new VisibilityResampler();
2683 : // Make the multi-threaded re-sampler and supply the basic
2684 : // re-sampler used in the worklet threads.
2685 : CountedPtr<VisibilityResamplerBase> mthVisResampler = new MultiThreadedVisibilityResampler(useDoublePrecGrid,
2686 : visResamplerCtor);
2687 : */
2688 : // Now make the FTMachine
2689 0 : if(facets_p>1) {
2690 : os << LogIO::NORMAL // Loglevel INFO
2691 : << "Multi-facet Fourier transforms will use specified common tangent point:"
2692 0 : << LogIO::POST;
2693 0 : os << LogIO::NORMAL << tangentPoint() << LogIO::POST; // Loglevel INFO
2694 : // ft_p = new rGridFT(cache_p / 2, tile_p, mthVisResampler, gridfunction_p, mLocation_p,
2695 : // phaseCenter_p, padding, false, useDoublePrecGrid);
2696 :
2697 0 : ft_p=new rGridFT(cache_p / 2, tile_p, gridfunction_p, mLocation_p,
2698 0 : phaseCenter_p, padding, false, useDoublePrecGrid);
2699 :
2700 : }
2701 : else {
2702 : os << LogIO::DEBUG1
2703 : << "Single facet Fourier transforms will use image center as tangent points"
2704 0 : << LogIO::POST;
2705 0 : ft_p = new rGridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p,
2706 0 : padding, false, useDoublePrecGrid);
2707 : }
2708 0 : AlwaysAssert(ft_p, AipsError);
2709 :
2710 0 : cft_p = new SimpleComponentFTMachine();
2711 0 : AlwaysAssert(cft_p, AipsError);
2712 :
2713 : }
2714 : //===============================================================
2715 : // A-Projection FTMachine code start here
2716 : //===============================================================
2717 37 : else if ((ftmachine_p == "awproject") || (ftmachine_p == "mawproject")){
2718 0 : if (wprojPlanes_p<=1)
2719 : {
2720 : os << LogIO::NORMAL
2721 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
2722 0 : << LogIO::POST;
2723 0 : os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
2724 : }
2725 0 : if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
2726 : {
2727 : os << LogIO::WARN
2728 : << "No. of w-planes set too low for W projection - recommend at least 128"
2729 0 : << LogIO::POST;
2730 0 : os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
2731 : }
2732 :
2733 0 : CountedPtr<ATerm> apertureFunction = createTelescopeATerm(*ms_p, aTermOn_p);
2734 0 : CountedPtr<PSTerm> psTerm = new PSTerm();
2735 0 : CountedPtr<WTerm> wTerm = new WTerm();
2736 :
2737 : //
2738 : // Selectively switch off CFTerms.
2739 : //
2740 0 : if (aTermOn_p == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
2741 0 : if (psTermOn_p == false) psTerm->setOpCode(CFTerms::NOOP);
2742 :
2743 : //
2744 : // Construct the CF object with appropriate CFTerms.
2745 : //
2746 0 : CountedPtr<ConvolutionFunction> awConvFunc;
2747 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP_p);
2748 : // if ((ftmachine_p=="mawproject") || (mTermOn_p))
2749 : // awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP_p);
2750 : // else
2751 0 : awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP_p);
2752 :
2753 : //
2754 : // Construct the appropriate re-sampler.
2755 : //
2756 0 : CountedPtr<VisibilityResamplerBase> visResampler = new AWVisResampler();
2757 : // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
2758 :
2759 : //
2760 : // Construct and initialize the CF cache object.
2761 : //
2762 0 : CountedPtr<CFCache> cfcache = new CFCache();
2763 0 : cfcache->setCacheDir(cfCacheDirName_p.data());
2764 0 : cfcache->initCache2();
2765 :
2766 : //
2767 : // Finally construct the FTMachine with the CFCache, ConvFunc and
2768 : // Re-sampler objects.
2769 : //
2770 0 : useDoublePrecGrid=(!singlePrec_p);
2771 0 : ft_p = new AWProjectWBFT(wprojPlanes_p, cache_p/2,
2772 : cfcache, awConvFunc,
2773 : // mthVisResampler,
2774 : visResampler,
2775 0 : /*true */doPointing, doPBCorr,
2776 0 : tile_p, computePAStep_p, pbLimit_p, true,conjBeams_p,
2777 0 : useDoublePrecGrid);
2778 :
2779 0 : ((AWProjectWBFT *)ft_p)->setObservatoryLocation(mLocation_p);
2780 : //
2781 : // Explicit type casting of ft_p does not look good. It does not
2782 : // pick up the setPAIncrement() method of PBWProjectFT without
2783 : // this
2784 : //
2785 : // os << LogIO::NORMAL << "Setting PA increment to " << parAngleInc_p.getValue("deg") << " deg" << endl;
2786 :
2787 : // ((AWProjectFT *)ft_p)->setPAIncrement(parAngleInc_p);
2788 0 : Quantity rotateOTF(rotPAStep_p,"deg");
2789 0 : ((AWProjectFT *)ft_p)->setPAIncrement(Quantity(computePAStep_p,"deg"),rotateOTF);
2790 :
2791 0 : AlwaysAssert(ft_p, AipsError);
2792 0 : cft_p = new SimpleComponentFTMachine();
2793 0 : AlwaysAssert(cft_p, AipsError);
2794 :
2795 0 : }
2796 : //===============================================================
2797 : // A-Projection FTMachine code end here
2798 : //===============================================================
2799 37 : else if(ftmachine_p=="SetJyGridFT"){
2800 0 : Vector<Double> freqs(0);
2801 0 : Vector<Double> scale(0);
2802 : os << LogIO::DEBUG1
2803 : << "SetJy Fourier transforms"
2804 0 : << LogIO::POST;
2805 0 : ft_p=new SetJyGridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p,
2806 0 : phaseCenter_p,
2807 0 : padding, false, useDoublePrecGrid, freqs, scale);
2808 0 : cft_p = new SimpleComponentFTMachine();
2809 :
2810 0 : }
2811 : else {
2812 : os << LogIO::NORMAL // Loglevel INFO
2813 : << "Performing interferometric gridding..."
2814 37 : << LogIO::POST;
2815 : os << LogIO::NORMAL1 // gridfunction_p is too cryptic for most users.
2816 37 : << "...with convolution function " << gridfunction_p << LogIO::POST;
2817 : // Now make the FTMachine
2818 37 : if(facets_p>1) {
2819 : os << LogIO::NORMAL // Loglevel INFO
2820 : << "Multi-facet Fourier transforms will use specified common tangent point:"
2821 0 : << LogIO::POST;
2822 0 : os << LogIO::NORMAL << tangentPoint() << LogIO::POST; // Loglevel INFO
2823 0 : ft_p = new GridFT(cache_p / 2, tile_p, gridfunction_p, mLocation_p,
2824 0 : phaseCenter_p, padding, false, useDoublePrecGrid);
2825 :
2826 : }
2827 : else {
2828 : os << LogIO::DEBUG1
2829 : << "Single facet Fourier transforms will use image center as tangent points"
2830 37 : << LogIO::POST;
2831 74 : ft_p = new GridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p,
2832 37 : padding, false, useDoublePrecGrid);
2833 :
2834 : }
2835 37 : AlwaysAssert(ft_p, AipsError);
2836 :
2837 37 : cft_p = new SimpleComponentFTMachine();
2838 37 : AlwaysAssert(cft_p, AipsError);
2839 :
2840 : }
2841 37 : ft_p->setnumthreads(numthreads_p);
2842 37 : cft_p->setnumthreads(numthreads_p);
2843 37 : ft_p->setSpw(dataspectralwindowids_p, freqFrameValid_p);
2844 37 : ft_p->setFreqInterpolation(freqInterpMethod_p);
2845 37 : if(doTrackSource_p){
2846 0 : ft_p->setMovingSource(trackDir_p);
2847 : }
2848 37 : ft_p->setSpwChanSelection(spwchansels_p);
2849 37 : ft_p->setSpwFreqSelection(mssFreqSel_p);
2850 :
2851 : /******* Start MTFT code ********/
2852 : // MultiTermFT is a container for an FTMachine of any type.
2853 : // It will apply Taylor-polynomial weights during gridding and degridding
2854 : // and will do multi-term grid-correction (normalizations).
2855 : // (1) ft_p already holds an FT of the correct type
2856 : // (2) If nterms>1, create a new MultiTermFT using ft_p, and reassign ft_p.
2857 : // Currently, Multi-Term applies only to wideband imaging.
2858 37 : if( ntaylor_p > 1 )
2859 : {
2860 : //cout << "UUU : Creating a Multi-Term FT machine containing " << ftmachine_p << endl;
2861 : FTMachine *tempftm;
2862 0 : if ( useNewMTFT_p == false )
2863 : {
2864 0 : tempftm = new MultiTermFT(ft_p, ft_p->name(), ntaylor_p, reffreq_p);
2865 : }
2866 : else
2867 : {
2868 0 : tempftm = new NewMultiTermFT(ft_p, ntaylor_p, reffreq_p);
2869 : }
2870 0 : ft_p = tempftm;
2871 : }
2872 : /******* End MTFT code ********/
2873 :
2874 37 : return true;
2875 37 : }
2876 :
2877 0 : String Imager::tangentPoint()
2878 : {
2879 0 : MVAngle mvRa=phaseCenter_p.getAngle().getValue()(0);
2880 0 : MVAngle mvDec=phaseCenter_p.getAngle().getValue()(1);
2881 0 : ostringstream oos;
2882 0 : oos << " ";
2883 0 : Int widthRA=20;
2884 0 : Int widthDec=20;
2885 0 : oos.setf(ios::left, ios::adjustfield);
2886 0 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
2887 0 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
2888 0 : oos << " "
2889 0 : << MDirection::showType(phaseCenter_p.getRefPtr()->getType());
2890 0 : return String(oos);
2891 0 : }
2892 :
2893 0 : Bool Imager::removeTable(const String& tablename) {
2894 :
2895 0 : LogIO os(LogOrigin("imager", "removeTable()", WHERE));
2896 :
2897 0 : if(Table::isReadable(tablename)) {
2898 0 : if (! Table::isWritable(tablename)) {
2899 : os << LogIO::SEVERE << "Table " << tablename
2900 0 : << " is not writable!: cannot alter it" << LogIO::POST;
2901 0 : return false;
2902 : }
2903 : else {
2904 0 : if (Table::isOpened(tablename)) {
2905 : os << LogIO::SEVERE << "Table " << tablename
2906 : << " is already open in the process. It needs to be closed first"
2907 0 : << LogIO::POST;
2908 0 : return false;
2909 : } else {
2910 0 : Table table(tablename, Table::Update);
2911 0 : if (table.isMultiUsed()) {
2912 : os << LogIO::SEVERE << "Table " << tablename
2913 : << " is already open in another process. It needs to be closed first"
2914 0 : << LogIO::POST;
2915 0 : return false;
2916 : } else {
2917 0 : Table table(tablename, Table::Delete);
2918 0 : }
2919 0 : }
2920 : }
2921 : }
2922 0 : return true;
2923 0 : }
2924 :
2925 0 : Bool Imager::updateSkyModel(const Vector<String>& model,
2926 : const String complist) {
2927 0 : LogIO os(LogOrigin("imager", "updateSkyModel()", WHERE));
2928 0 : if(redoSkyModel_p)
2929 0 : throw(AipsError("Programming error: update skymodel is called without a valid skymodel"));
2930 0 : Bool coordMatch=true;
2931 0 : for (Int thismodel=0;thismodel<Int(model.nelements());++thismodel) {
2932 0 : CoordinateSystem cs=(sm_p->image(thismodel)).coordinates();
2933 0 : coordMatch= coordMatch || checkCoord(cs, model(thismodel));
2934 : ///return false if any fails anyways
2935 0 : if(!coordMatch)
2936 0 : return false;
2937 0 : if(model(thismodel)=="") {
2938 : os << LogIO::SEVERE << "Need a name for model "
2939 0 : << model << LogIO::POST;
2940 0 : return false;
2941 : }
2942 : else {
2943 0 : if(!Table::isReadable(model(thismodel))) {
2944 : os << LogIO::SEVERE << model(thismodel) << "is unreadable"
2945 0 : << model << LogIO::POST;
2946 0 : return false;
2947 : }
2948 : }
2949 0 : images_p[thismodel]=0;
2950 0 : images_p[thismodel]=new PagedImage<Float>(model(thismodel));
2951 0 : AlwaysAssert(!images_p[thismodel].null(), AipsError);
2952 0 : sm_p->updatemodel(thismodel, *images_p[thismodel]);
2953 0 : }
2954 0 : if((complist !="") && Table::isReadable(complist)){
2955 0 : ComponentList cl(Path(complist), true);
2956 0 : sm_p->updatemodel(cl);
2957 0 : }
2958 0 : return true;
2959 0 : }
2960 :
2961 0 : Bool Imager::createSkyEquation(const String complist)
2962 : {
2963 0 : Vector<String> image;
2964 0 : Vector<String> mask;
2965 0 : Vector<String> fluxMask;
2966 0 : Vector<Bool> fixed;
2967 0 : return createSkyEquation(image, fixed, mask, fluxMask, complist);
2968 0 : }
2969 :
2970 37 : Bool Imager::createSkyEquation(const Vector<String>& image,
2971 : const String complist)
2972 : {
2973 37 : Vector<Bool> fixed(image.nelements()); fixed=false;
2974 37 : Vector<String> mask(image.nelements()); mask="";
2975 37 : Vector<String> fluxMask(image.nelements()); fluxMask="";
2976 74 : return createSkyEquation(image, fixed, mask, fluxMask, complist);
2977 37 : }
2978 :
2979 0 : Bool Imager::createSkyEquation(const Vector<String>& image, const Vector<Bool>& fixed,
2980 : const String complist)
2981 : {
2982 0 : Vector<String> mask(image.nelements()); mask="";
2983 0 : Vector<String> fluxMask(image.nelements()); fluxMask="";
2984 0 : return createSkyEquation(image, fixed, mask, fluxMask, complist);
2985 0 : }
2986 :
2987 0 : Bool Imager::createSkyEquation(const Vector<String>& image, const Vector<Bool>& fixed,
2988 : const Vector<String>& mask,
2989 : const String complist)
2990 : {
2991 0 : Vector<String> fluxMask(image.nelements()); fluxMask="";
2992 0 : return createSkyEquation(image, fixed, mask, fluxMask, complist);
2993 0 : }
2994 :
2995 37 : Bool Imager::createSkyEquation(const Vector<String>& image,
2996 : const Vector<Bool>& fixed,
2997 : const Vector<String>& mask,
2998 : const Vector<String>& fluxMask,
2999 : const String complist)
3000 : {
3001 :
3002 37 : if(!valid()) return false;
3003 :
3004 74 : LogIO os(LogOrigin("imager", "createSkyEquation()", WHERE));
3005 :
3006 : // If there is no sky model, we'll make one:
3007 :
3008 37 : if(sm_p==0) {
3009 37 : if((facets_p >1)){
3010 : // Support serial and parallel specializations
3011 0 : setWFCleanImageSkyModel();
3012 : }
3013 37 : else if(ntaylor_p>1){
3014 : // Init the msmfs sky-model so that Taylor weights are triggered in CubeSkyEqn
3015 0 : sm_p = new WBCleanImageSkyModel(ntaylor_p, 1 ,reffreq_p);
3016 : }
3017 : else {
3018 37 : sm_p = new CleanImageSkyModel();
3019 : }
3020 : }
3021 37 : AlwaysAssert(sm_p, AipsError);
3022 : //Now lets tell it how to use memory and templattices
3023 37 : sm_p->setMemoryUse(avoidTempLatt_p);
3024 37 : if(imageTileVol_p > 1000)
3025 0 : sm_p->setTileVol(imageTileVol_p);
3026 :
3027 : // Read data pol rep from the MS. This is also done in imagecoordinates().
3028 : // This code segment is repeated here because imagecoordinates
3029 : // is not always called before SkyModel is made (im.ft).
3030 : // Note : This should go into a function.
3031 : {
3032 37 : MSColumns msc(*ms_p);
3033 37 : Vector<String> polType=msc.feed().polarizationType()(0);
3034 111 : if (polType(0)!="X" && polType(0)!="Y" &&
3035 111 : polType(0)!="R" && polType(0)!="L") {
3036 : os << LogIO::WARN << "Unknown stokes types in feed table: ["
3037 0 : << polType(0) << ", " << polType(1) << "]" << endl
3038 0 : << "Results open to question!" << LogIO::POST;
3039 : }
3040 :
3041 37 : if (polType(0)=="X" || polType(0)=="Y") {
3042 0 : polRep_p=StokesImageUtil::LINEAR;
3043 : os << LogIO::DEBUG1
3044 0 : << "Preferred polarization representation is linear" << LogIO::POST;
3045 : }
3046 : else {
3047 37 : polRep_p=StokesImageUtil::CIRCULAR;
3048 : os << LogIO::DEBUG1
3049 37 : << "Preferred polarization representation is circular" << LogIO::POST;
3050 : }
3051 37 : }// end of reading-in polRep_p from the MS
3052 :
3053 : // Send the data correlation type to the SkyModel
3054 : // This needs to be done before the first call to 'ImageSkyModel::cImage()'
3055 37 : os << LogIO::DEBUG1 << "Data PolRep in Imager2.cc::createSkyEquation : " << polRep_p << LogIO::POST;
3056 37 : sm_p->setDataPolFrame(polRep_p);
3057 :
3058 :
3059 : // Add the componentlist
3060 : // cerr << "COMPLIST " << complist << endl;
3061 37 : if(complist!="") {
3062 37 : if(!Table::isReadable(complist)) {
3063 : os << LogIO::SEVERE << "ComponentList " << complist
3064 0 : << " not readable" << LogIO::POST;
3065 0 : return false;
3066 : }
3067 37 : if(componentList_p){
3068 0 : delete componentList_p;
3069 0 : componentList_p=0;
3070 : }
3071 37 : componentList_p=new ComponentList(complist, true);
3072 37 : if(componentList_p==0) {
3073 : os << LogIO::SEVERE << "Cannot create ComponentList from " << complist
3074 0 : << LogIO::POST;
3075 0 : return false;
3076 : }
3077 37 : if(!sm_p->add(*componentList_p)) {
3078 : os << LogIO::SEVERE << "Cannot add ComponentList " << complist
3079 0 : << " to SkyModel" << LogIO::POST;
3080 0 : return false;
3081 : }
3082 : os << LogIO::NORMAL // Loglevel INFO
3083 37 : << "Processing after subtracting componentlist " << complist << LogIO::POST;
3084 : }
3085 : else {
3086 0 : delete componentList_p;
3087 0 : componentList_p=0;
3088 : }
3089 :
3090 : // Make image with the required shape and coordinates only if
3091 : // they don't exist yet
3092 37 : nmodels_p=image.nelements();
3093 :
3094 : // Remove degenerate case (due to user interface?)
3095 37 : if((nmodels_p==1)&&(image(0)=="")) {
3096 0 : nmodels_p=0;
3097 : }
3098 37 : if(nmodels_p>0) {
3099 0 : images_p.resize(nmodels_p);
3100 0 : masks_p.resize(nmodels_p);
3101 0 : fluxMasks_p.resize(nmodels_p);
3102 0 : residuals_p.resize(nmodels_p);
3103 0 : for (Int model=0;model<Int(nmodels_p);++model) {
3104 0 : if(image(model)=="") {
3105 : os << LogIO::SEVERE << "Need a name for model "
3106 0 : << model << LogIO::POST;
3107 0 : return false;
3108 : }
3109 : else {
3110 0 : if(!Table::isReadable(image(model))) {
3111 0 : if(!assertDefinedImageParameters()) return false;
3112 0 : make(image(model));
3113 : }
3114 : }
3115 0 : images_p[model]=0;
3116 0 : images_p[model]=new PagedImage<Float>(image(model));
3117 0 : AlwaysAssert(!images_p[model].null(), AipsError);
3118 :
3119 : //Determining the number of XFR
3120 0 : Int numOfXFR=nmodels_p+1;
3121 0 : if(datafieldids_p.nelements() >0)
3122 0 : numOfXFR=datafieldids_p.nelements()*nmodels_p + 1;
3123 0 : if(squintType_p != BeamSquint::NONE){
3124 0 : if(parAngleInc_p.getValue("deg") >0 ){
3125 0 : numOfXFR= numOfXFR* Int(360/parAngleInc_p.getValue("deg"));
3126 : }
3127 : else{
3128 0 : numOfXFR= numOfXFR*10;
3129 : }
3130 : }
3131 0 : if((sm_p->add(*images_p[model], numOfXFR))!=model) {
3132 0 : os << LogIO::SEVERE << "Error adding model " << model << LogIO::POST;
3133 0 : return false;
3134 : }
3135 :
3136 0 : fluxMasks_p[model]=0;
3137 0 : if(fluxMask(model)!=""&&Table::isReadable(fluxMask(model))) {
3138 0 : fluxMasks_p[model]=new PagedImage<Float>(fluxMask(model));
3139 0 : AlwaysAssert(!fluxMasks_p[model].null(), AipsError);
3140 0 : if(!sm_p->addFluxMask(model, *fluxMasks_p[model])) {
3141 : os << LogIO::SEVERE << "Error adding flux mask " << model
3142 0 : << " : " << fluxMask(model) << LogIO::POST;
3143 0 : return false;
3144 : }
3145 : }
3146 0 : residuals_p[model]=0;
3147 : }
3148 0 : addMasksToSkyEquation(mask, fixed);
3149 : }
3150 :
3151 : // Always need a VisSet and an FTMachine
3152 37 : if (!ft_p)
3153 37 : createFTMachine();
3154 :
3155 : // Now set up the SkyEquation
3156 37 : AlwaysAssert(sm_p, AipsError);
3157 : // AlwaysAssert(vs_p, AipsError);
3158 37 : AlwaysAssert(rvi_p, AipsError);
3159 37 : AlwaysAssert(ft_p, AipsError);
3160 37 : AlwaysAssert(cft_p, AipsError);
3161 :
3162 : // Setup the sky equation
3163 37 : setSkyEquation();
3164 :
3165 : // If primary-beams are needed, force the fluxScale images held by
3166 : // the SkyModel classes to be allocated/initialized.
3167 37 : if(doVP_p){
3168 0 : if( (sm_p->numberOfModels() > 0) && (ft_p->name() != "MosaicFT"))
3169 0 : sm_p->mandateFluxScale(0);
3170 : }
3171 :
3172 :
3173 : /* // Commented out by URV (4 Apr 2012) : Redundant Code.
3174 : // This block determines which SkyEquation is to be used.
3175 : // We are using a mf* algorithm and there is more than one image
3176 :
3177 : if (doMultiFields_p && multiFields_p) {
3178 : // Mosaicing case
3179 : if(doVP_p){
3180 : //bypassing the minimum size FT stuff as its broken till its fixed
3181 : //se_p=new MosaicSkyEquation(*sm_p, *vs_p, *ft_p, *cft_p);
3182 :
3183 : setSkyEquation();
3184 : if(ft_p->name() != "MosaicFT")
3185 : sm_p->mandateFluxScale(0);
3186 : os << LogIO::NORMAL // Loglevel INFO
3187 : << "Mosaicing multiple fields with simple sky equation" << LogIO::POST;
3188 : }
3189 : // mosaicing with no vp correction
3190 : else{
3191 : setSkyEquation();
3192 : os << LogIO::NORMAL // Loglevel INFO
3193 : << "Processing multiple fields with simple sky equation" << LogIO::POST;
3194 : os << LogIO::WARN
3195 : << "Voltage Pattern is not set: will not correct for primary beam"
3196 : << LogIO::POST;
3197 : doMultiFields_p=false;
3198 : }
3199 : }
3200 : // We are not using an mf* algorithm or there is only one image
3201 : else {
3202 : // Support serial and parallel specializations
3203 : if((facets_p >1)){
3204 : setSkyEquation();
3205 : //se_p=new SkyEquation(*sm_p, *vs_p, *ft_p, *cft_p);
3206 : os << LogIO::NORMAL // Loglevel INFO
3207 : << "Processing multiple facets with simple sky equation" << LogIO::POST;
3208 : }
3209 : // Mosaicing
3210 : else if(doVP_p) {
3211 : //Bypassing the mosaicskyequation to the slow version for now.
3212 : // se_p=new MosaicSkyEquation(*sm_p, *vs_p, *ft_p, *cft_p);
3213 :
3214 : setSkyEquation();
3215 : if(ft_p->name() != "MosaicFT")
3216 : sm_p->mandateFluxScale(0);
3217 : os << LogIO::NORMAL // Loglevel PROGRESS
3218 : << "Mosaicing single field with simple sky equation" << LogIO::POST;
3219 : }
3220 : // Default
3221 : else {
3222 : setSkyEquation();
3223 : os << LogIO::DEBUG1
3224 : << "Processing single field with simple sky equation" << LogIO::POST;
3225 : }
3226 : }
3227 :
3228 : */
3229 :
3230 :
3231 : //os.localSink().flush();
3232 : //For now force none sault weighting with mosaic ft machine
3233 :
3234 37 : String scaleType=scaleType_p;
3235 37 : if(ft_p->name()=="MosaicFT")
3236 0 : scaleType="NONE";
3237 :
3238 37 : se_p->setImagePlaneWeighting(scaleType, minPB_p, constPB_p);
3239 37 : se_p->doFlatNoise(flatnoise_p);
3240 :
3241 37 : AlwaysAssert(se_p, AipsError);
3242 :
3243 : // Now add any SkyJones that are needed
3244 37 : if(doVP_p && (ft_p->name()!="MosaicFT")) {
3245 0 : MSColumns msc(*ms_p);
3246 0 : if (doDefaultVP_p) {
3247 0 : vp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p, skyPosThreshold_p);
3248 : } else {
3249 0 : Table vpTable( vpTableStr_p );
3250 0 : vp_p=new VPSkyJones(msc, vpTable, parAngleInc_p, squintType_p, skyPosThreshold_p);
3251 0 : }
3252 0 : se_p->setSkyJones(*vp_p);
3253 0 : }
3254 : else {
3255 37 : vp_p=0;
3256 : }
3257 37 : return true;
3258 37 : }
3259 :
3260 0 : Bool Imager::addResiduals(const Vector<String>& imageNames) {
3261 0 : Bool retval=true;
3262 0 : residuals_p.resize(imageNames.nelements(), true, false);
3263 0 : for (Int thismodel=0;thismodel<Int(imageNames.nelements());++thismodel) {
3264 0 : if(imageNames(thismodel)!="") {
3265 0 : residuals_p[thismodel]=0;
3266 0 : if(Table::isWritable(imageNames(thismodel))){
3267 0 : residuals_p[thismodel]=new PagedImage<Float>(imageNames(thismodel));
3268 0 : if(!(residuals_p[thismodel]->shape()).isEqual(images_p[thismodel]->shape())){
3269 0 : residuals_p[thismodel]=0;
3270 0 : removeTable(imageNames(thismodel));
3271 : }
3272 : }
3273 0 : if(residuals_p[thismodel].null()){
3274 0 : if(Table::isReadable(imageNames(thismodel)))
3275 0 : removeTable(imageNames(thismodel));
3276 0 : residuals_p[thismodel]=
3277 0 : new PagedImage<Float> (TiledShape(images_p[thismodel]->shape(),
3278 0 : images_p[thismodel]->niceCursorShape()),
3279 0 : images_p[thismodel]->coordinates(),
3280 0 : imageNames(thismodel));
3281 0 : AlwaysAssert(!residuals_p[thismodel].null(), AipsError);
3282 0 : residuals_p[thismodel]->setUnits(Unit("Jy/beam"));
3283 : }
3284 0 : if(residuals_p[thismodel].null())
3285 0 : retval=false;
3286 : }
3287 : else{
3288 0 : retval=false;
3289 : }
3290 : }
3291 0 : return retval;
3292 : }
3293 : // Tell the sky model to use the specified images as the residuals
3294 0 : Bool Imager::addResidualsToSkyEquation(const Vector<String>& imageNames) {
3295 :
3296 0 : addResiduals(imageNames);
3297 0 : for (Int thismodel=0;thismodel<Int(imageNames.nelements());++thismodel) {
3298 0 : if(imageNames(thismodel)!="")
3299 0 : sm_p->addResidual(thismodel, *residuals_p[thismodel]);
3300 : }
3301 0 : return true;
3302 : }
3303 :
3304 112 : void Imager::destroySkyEquation()
3305 : {
3306 112 : if(se_p) delete se_p; se_p=0;
3307 112 : if(sm_p) delete sm_p; sm_p=0;
3308 112 : if(vp_p) delete vp_p; vp_p=0;
3309 112 : if(gvp_p) delete gvp_p; gvp_p=0;
3310 :
3311 112 : if(componentList_p) delete componentList_p; componentList_p=0;
3312 :
3313 :
3314 112 : for (Int model=0;model<Int(nmodels_p); ++model) {
3315 : //As these are CountedPtrs....just assigning them to NULL
3316 : //get the objects out of context
3317 0 : if(!images_p[model].null())
3318 0 : images_p[model]=0;
3319 :
3320 0 : if(!masks_p[model].null())
3321 0 : masks_p[model]=0;
3322 :
3323 0 : if(!fluxMasks_p[model].null())
3324 0 : fluxMasks_p[model]=0;
3325 :
3326 0 : if(!residuals_p[model].null())
3327 0 : residuals_p[model]=0;
3328 : }
3329 :
3330 112 : redoSkyModel_p=true;
3331 112 : }
3332 :
3333 0 : Bool Imager::assertDefinedImageParameters() const
3334 : {
3335 0 : LogIO os(LogOrigin("imager", "if(!assertDefinedImageParameters()", WHERE));
3336 0 : if(!setimaged_p) {
3337 : os << LogIO::SEVERE << "Image parameters not yet set: use im.defineimage."
3338 0 : << LogIO::POST;
3339 0 : return false;
3340 : }
3341 0 : return true;
3342 0 : }
3343 :
3344 127 : Bool Imager::valid() const {
3345 254 : LogIO os(LogOrigin("imager", "if(!valid()) return false", WHERE));
3346 127 : if(ms_p.null()) {
3347 : os << LogIO::SEVERE << "Program logic error: MeasurementSet pointer ms_p not yet set"
3348 0 : << LogIO::POST;
3349 0 : return false;
3350 : }
3351 127 : if(mssel_p.null()) {
3352 : os << LogIO::SEVERE << "Program logic error: MeasurementSet pointer mssel_p not yet set"
3353 0 : << LogIO::POST;
3354 0 : return false;
3355 : }
3356 127 : if(!rvi_p) {
3357 : os << LogIO::SEVERE << "Program logic error: VisibilityIterator not yet set"
3358 0 : << LogIO::POST;
3359 0 : return false;
3360 : }
3361 127 : return true;
3362 127 : }
3363 :
3364 :
3365 0 : Bool Imager::addMasksToSkyEquation(const Vector<String>& mask, const Vector<Bool>& fixed){
3366 0 : LogIO os(LogOrigin("imager", "addMasksToSkyEquation()", WHERE));
3367 :
3368 0 : for(Int model=0 ;model < nmodels_p; ++model){
3369 :
3370 :
3371 0 : if((Int(fixed.nelements())>model) && fixed(model)) {
3372 : os << LogIO::NORMAL // Loglevel INFO
3373 0 : << "Model " << model << " will be held fixed" << LogIO::POST;
3374 0 : sm_p->fix(model);
3375 : }
3376 : /*
3377 : if(!(masks_p[model].null())) delete masks_p[model];
3378 : masks_p[model]=0;
3379 : */
3380 0 : if(mask(model)!=""&&Table::isReadable(mask(model))) {
3381 0 : masks_p[model]=new PagedImage<Float>(mask(model));
3382 0 : AlwaysAssert(!masks_p[model].null(), AipsError);
3383 0 : if(!sm_p->addMask(model, *masks_p[model])) {
3384 : os << LogIO::SEVERE << "Error adding mask " << model
3385 0 : << " : " << mask(model) << LogIO::POST;
3386 0 : return false;
3387 : }
3388 : }
3389 : }
3390 0 : return true;
3391 0 : }
3392 :
3393 : void
3394 112 : Imager::openSubTable (const Table & otherTable, Table & table, const TableLock & tableLock)
3395 : {
3396 112 : if (otherTable.isNull()){
3397 :
3398 : // otherTable does not exist so leave things be
3399 :
3400 : }
3401 88 : else if (otherTable.tableType() == Table::Memory){
3402 :
3403 0 : table = otherTable;
3404 :
3405 : }
3406 : else{
3407 :
3408 : // Reopen (potentially) the subtable with the desired locking
3409 :
3410 88 : table = Table (otherTable.tableName(), tableLock);
3411 : }
3412 112 : }
3413 :
3414 : Bool
3415 7 : Imager::openSubTables()
3416 : {
3417 : // These variables will already have copied in the Tables from
3418 : // the MS specified in open. If they are not memory resident
3419 : // subtables then replace them with table objects having the
3420 : // UserNoReadLocking attribute.
3421 :
3422 7 : TableLock tableLock (TableLock::UserNoReadLocking);
3423 :
3424 7 : openSubTable (ms_p->antenna(), antab_p, tableLock);
3425 7 : openSubTable (ms_p->dataDescription (), datadesctab_p, tableLock);
3426 7 : openSubTable (ms_p->doppler(), dopplertab_p, tableLock);
3427 7 : openSubTable (ms_p->feed(), feedtab_p, tableLock);
3428 7 : openSubTable (ms_p->field(), fieldtab_p, tableLock);
3429 7 : openSubTable (ms_p->flagCmd(), flagcmdtab_p, tableLock);
3430 7 : openSubTable (ms_p->freqOffset(), freqoffsettab_p, tableLock);
3431 7 : openSubTable (ms_p->observation(), obstab_p, tableLock);
3432 7 : openSubTable (ms_p->pointing(), pointingtab_p, tableLock);
3433 7 : openSubTable (ms_p->polarization(), poltab_p, tableLock);
3434 7 : openSubTable (ms_p->processor(), proctab_p, tableLock);
3435 7 : openSubTable (ms_p->source(), sourcetab_p, tableLock);
3436 7 : openSubTable (ms_p->spectralWindow(), spwtab_p, tableLock);
3437 7 : openSubTable (ms_p->state(), statetab_p, tableLock);
3438 7 : openSubTable (ms_p->sysCal(), syscaltab_p, tableLock);
3439 7 : openSubTable (ms_p->weather(), weathertab_p, tableLock);
3440 :
3441 : // Handle the history table
3442 :
3443 7 : if(ms_p->isWritable()){
3444 :
3445 7 : if(!(Table::isReadable(ms_p->historyTableName()))){
3446 :
3447 : // setup a new table in case its not there
3448 0 : TableRecord &kws = ms_p->rwKeywordSet();
3449 0 : SetupNewTable historySetup(ms_p->historyTableName(),
3450 0 : MSHistory::requiredTableDesc(),Table::New);
3451 0 : kws.defineTable(MS::keywordName(MS::HISTORY), Table(historySetup));
3452 :
3453 0 : }
3454 21 : historytab_p=Table(ms_p->historyTableName(),
3455 14 : TableLock(TableLock::UserNoReadLocking), Table::Update);
3456 :
3457 7 : hist_p= new MSHistoryHandler(*ms_p, "imager");
3458 : }
3459 :
3460 7 : return true;
3461 :
3462 : }
3463 :
3464 97 : Bool Imager::lock(){
3465 :
3466 : Bool ok;
3467 97 : ok=true;
3468 97 : if(lockCounter_p == 0){
3469 :
3470 14 : ok= ok && (ms_p->lock());
3471 14 : ok= ok && antab_p.lock(false);
3472 14 : ok= ok && datadesctab_p.lock(false);
3473 14 : ok= ok && feedtab_p.lock(false);
3474 14 : ok= ok && fieldtab_p.lock(false);
3475 14 : ok= ok && obstab_p.lock(false);
3476 14 : ok= ok && poltab_p.lock(false);
3477 14 : ok= ok && proctab_p.lock(false);
3478 14 : ok= ok && spwtab_p.lock(false);
3479 14 : ok= ok && statetab_p.lock(false);
3480 14 : if(!dopplertab_p.isNull())
3481 0 : ok= ok && dopplertab_p.lock(false);
3482 14 : if(!flagcmdtab_p.isNull())
3483 14 : ok= ok && flagcmdtab_p.lock(false);
3484 14 : if(!freqoffsettab_p.isNull())
3485 0 : ok= ok && freqoffsettab_p.lock(false);
3486 14 : if(!historytab_p.isNull())
3487 14 : ok= ok && historytab_p.lock(false);
3488 14 : if(!pointingtab_p.isNull())
3489 14 : ok= ok && pointingtab_p.lock(false);
3490 14 : if(!sourcetab_p.isNull())
3491 14 : ok= ok && sourcetab_p.lock(false);
3492 14 : if(!syscaltab_p.isNull())
3493 4 : ok= ok && syscaltab_p.lock(false);
3494 14 : if(!weathertab_p.isNull())
3495 4 : ok= ok && weathertab_p.lock(false);
3496 :
3497 : }
3498 97 : ++lockCounter_p;
3499 :
3500 97 : return ok ;
3501 : }
3502 :
3503 126 : Bool Imager::unlock(){
3504 :
3505 126 : if(lockCounter_p==1){
3506 14 : ms_p->unlock();
3507 14 : antab_p.unlock();
3508 14 : datadesctab_p.unlock();
3509 14 : feedtab_p.unlock();
3510 14 : fieldtab_p.unlock();
3511 14 : obstab_p.unlock();
3512 14 : poltab_p.unlock();
3513 14 : proctab_p.unlock();
3514 14 : spwtab_p.unlock();
3515 14 : statetab_p.unlock();
3516 14 : if(!dopplertab_p.isNull())
3517 0 : dopplertab_p.unlock();
3518 14 : if(!flagcmdtab_p.isNull())
3519 14 : flagcmdtab_p.unlock();
3520 14 : if(!freqoffsettab_p.isNull())
3521 0 : freqoffsettab_p.unlock();
3522 14 : if(!historytab_p.isNull())
3523 14 : historytab_p.unlock();
3524 14 : if(!pointingtab_p.isNull())
3525 14 : pointingtab_p.unlock();
3526 14 : if(!sourcetab_p.isNull())
3527 14 : sourcetab_p.unlock();
3528 14 : if(!syscaltab_p.isNull())
3529 4 : syscaltab_p.unlock();
3530 14 : if(!weathertab_p.isNull())
3531 4 : weathertab_p.unlock();
3532 : }
3533 126 : for (Int thismodel=0;thismodel<Int(images_p.nelements());++thismodel) {
3534 0 : if ((images_p.nelements() > uInt(thismodel)) && (!images_p[thismodel].null())) {
3535 0 : images_p[thismodel]->table().relinquishAutoLocks(true);
3536 0 : images_p[thismodel]->table().unlock();
3537 : }
3538 0 : if ((residuals_p.nelements()> uInt(thismodel)) && (!residuals_p[thismodel].null())) {
3539 0 : residuals_p[thismodel]->table().relinquishAutoLocks(true);
3540 0 : residuals_p[thismodel]->table().unlock();
3541 : }
3542 0 : if ((masks_p.nelements()> uInt(thismodel)) && (!masks_p[thismodel].null())) {
3543 0 : masks_p[thismodel]->table().relinquishAutoLocks(true);
3544 0 : masks_p[thismodel]->table().unlock();
3545 : }
3546 : }
3547 126 : if(lockCounter_p > 0 )
3548 97 : --lockCounter_p;
3549 126 : return true ;
3550 : }
3551 :
3552 46 : Bool Imager::selectDataChannel(Vector<Int>& spectralwindowids,
3553 : String& dataMode,
3554 : Vector<Int>& dataNchan,
3555 : Vector<Int>& dataStart, Vector<Int>& dataStep,
3556 : MRadialVelocity& /*mDataStart*/,
3557 : MRadialVelocity& /*mDataStep*/){
3558 :
3559 :
3560 :
3561 92 : LogIO os(LogOrigin("Imager", "selectDataChannel()", WHERE));
3562 :
3563 :
3564 46 : if(dataMode=="channel") {
3565 46 : if (dataNchan.nelements() != spectralwindowids.nelements()){
3566 0 : if(dataNchan.nelements()==1){
3567 0 : dataNchan.resize(spectralwindowids.nelements(), true);
3568 0 : for(uInt k=1; k < spectralwindowids.nelements(); ++k){
3569 0 : dataNchan[k]=dataNchan[0];
3570 : }
3571 : }
3572 : else{
3573 : os << LogIO::SEVERE
3574 : << "Vector of nchan has to be of size 1 or be of the same shape as spw "
3575 0 : << LogIO::POST;
3576 0 : return false;
3577 : }
3578 : }
3579 46 : if (dataStart.nelements() != spectralwindowids.nelements()){
3580 0 : if(dataStart.nelements()==1){
3581 0 : dataStart.resize(spectralwindowids.nelements(), true);
3582 0 : for(uInt k=1; k < spectralwindowids.nelements(); ++k){
3583 0 : dataStart[k]=dataStart[0];
3584 : }
3585 : }
3586 : else{
3587 : os << LogIO::SEVERE
3588 : << "Vector of start has to be of size 1 or be of the same shape as spw "
3589 0 : << LogIO::POST;
3590 0 : return false;
3591 : }
3592 : }
3593 46 : if (dataStep.nelements() != spectralwindowids.nelements()){
3594 0 : if(dataStep.nelements()==1){
3595 0 : dataStep.resize(spectralwindowids.nelements(), true);
3596 0 : for(uInt k=1; k < spectralwindowids.nelements(); ++k){
3597 0 : dataStep[k]=dataStep[0];
3598 : }
3599 : }
3600 : else{
3601 : os << LogIO::SEVERE
3602 : << "Vector of step has to be of size 1 or be of the same shape as spw "
3603 0 : << LogIO::POST;
3604 0 : return false;
3605 : }
3606 : }
3607 :
3608 46 : if(spectralwindowids.nelements()>0) {
3609 46 : Int nch=0;
3610 152 : for(uInt i=0;i<spectralwindowids.nelements();++i) {
3611 106 : Int spwid=spectralwindowids(i);
3612 106 : Int numberChan=rvi_p->msColumns().spectralWindow().numChan()(spwid);
3613 106 : if(dataStart[i]<0) {
3614 : os << LogIO::SEVERE << "Illegal start pixel = "
3615 0 : << dataStart[i] << " for spw " << spwid
3616 0 : << LogIO::POST;
3617 0 : return false;
3618 : }
3619 :
3620 106 : if(dataNchan[i]<=0){
3621 37 : if(dataStep[i] <= 0)
3622 0 : dataStep[i]=1;
3623 37 : nch=(numberChan-dataStart[i])/Int(dataStep[i])+1;
3624 : }
3625 69 : else nch = dataNchan[i];
3626 143 : while((nch*dataStep[i]+dataStart[i]) > numberChan){
3627 37 : --nch;
3628 : }
3629 106 : Int end = Int(dataStart[i]) + Int(nch-1) * Int(dataStep[i]);
3630 106 : if(end < 0 || end > (numberChan)-1) {
3631 0 : os << LogIO::SEVERE << "Illegal step pixel = " << dataStep[i]
3632 : << " for spw " << spwid
3633 : << "\n end channel " << end
3634 0 : << " is out of range " << dataStart[i] << " to "
3635 : << (numberChan-1)
3636 0 : << LogIO::POST;
3637 0 : return false;
3638 : }
3639 :
3640 : os << LogIO::DEBUG1 // Too contentious for DEBUG1
3641 106 : << "Selecting within ";
3642 106 : if(nch > 1)
3643 : os << nch << " channels, starting at "
3644 106 : << dataStart[i] << ", stepped by " << dataStep[i] << ",";
3645 : else
3646 0 : os << "channel " << dataStart[i];
3647 106 : os << " for spw " << spwid << LogIO::POST;
3648 :
3649 : ///////////This is totally funked ...does not respect the spw selection
3650 : //whatever you do the the ngroups is always all the spw in the ms !!!
3651 : //vi.allSelectedSpectralWindows gets borked because of that
3652 : //rvi_p->selectChannel(1, Int(dataStart[i]), Int(nch),
3653 : // Int(dataStep[i]), spwid);
3654 106 : dataNchan[i]=nch;
3655 : }
3656 : /////Temporary replacement via the multims one
3657 46 : Block<Vector<Int> > blspw(1);
3658 46 : Block<Vector<Int> > blngr(1);
3659 46 : Block<Vector<Int> > blstart(1);
3660 46 : Block<Vector<Int> > blwid(1);
3661 46 : Block<Vector<Int> > blinr(1);
3662 46 : blspw[0]=spectralwindowids;
3663 46 : blngr[0]=Vector<Int>(spectralwindowids.nelements(),1);
3664 46 : blstart[0]=dataStart;
3665 46 : blwid=dataNchan;
3666 46 : blinr[0]=dataStep;
3667 46 : rvi_p->selectChannel(blngr, blstart, blwid,
3668 : blinr, blspw);
3669 : ////////////////////////
3670 :
3671 46 : } else {
3672 0 : VisBufferAutoPtr vb (rvi_p);
3673 0 : rvi_p->originChunks ();
3674 0 : Int numberChan=vb->msColumns().spectralWindow().numChan()(0);
3675 :
3676 0 : if(dataNchan[0]<=0){
3677 0 : if(dataStep[0] <=0)
3678 0 : dataStep[0]=1;
3679 0 : dataNchan[0]=(numberChan-dataStart[0])/Int(dataStep[0])+1;
3680 :
3681 : }
3682 0 : while((dataNchan[0]*dataStep[0]+dataStart[0]) > numberChan)
3683 0 : --dataNchan[0];
3684 :
3685 0 : Int end = Int(dataStart[0]) + Int(dataNchan[0]-1)
3686 0 : * Int(dataStep[0]);
3687 0 : if(end < 0 || end > (numberChan)-1) {
3688 0 : os << LogIO::SEVERE << "Illegal step pixel = " << dataStep[0]
3689 : << "\n end channel " << end << " is out of range 1 to "
3690 : << (numberChan-1)
3691 0 : << LogIO::POST;
3692 0 : return false;
3693 : }
3694 0 : os << LogIO::DEBUG1 << "Selecting within "<< dataNchan[0] // Loglevel INFO
3695 : << " channels, starting at visibility channel "
3696 0 : << dataStart[0] << " stepped by "
3697 0 : << dataStep[0] << LogIO::POST;
3698 0 : }
3699 : }
3700 :
3701 :
3702 46 : return true;
3703 :
3704 46 : }
3705 :
3706 :
3707 0 : Bool Imager::checkCoord(const CoordinateSystem& coordsys,
3708 : const String& imageName){
3709 :
3710 0 : PagedImage<Float> image(imageName);
3711 0 : CoordinateSystem imageCoord= image.coordinates();
3712 0 : Vector<Int> imageShape= image.shape().asVector();
3713 :
3714 0 : if(imageShape.nelements() > 3){
3715 0 : if(imageShape(3) != imageNchan_p)
3716 0 : return false;
3717 : }
3718 : else{
3719 0 : if(imageNchan_p >1)
3720 0 : return false;
3721 : }
3722 :
3723 0 : if(imageShape.nelements() > 2){
3724 0 : if(imageShape(2) != npol_p)
3725 0 : return false;
3726 : }
3727 : else{
3728 0 : if(npol_p > 1)
3729 0 : return false;
3730 : }
3731 0 : if(imageShape(0) != nx_p)
3732 0 : return false;
3733 0 : if(imageShape(1) != ny_p)
3734 0 : return false;
3735 :
3736 :
3737 :
3738 0 : if(!imageCoord.near(coordsys)){
3739 0 : return false;
3740 : }
3741 :
3742 : /*
3743 : DirectionCoordinate dir1(coordsys.directionCoordinate(0));
3744 : DirectionCoordinate dir2(imageCoord.directionCoordinate(0));
3745 : if(dir1.increment()(0) != dir2.increment()(0))
3746 : return false;
3747 : if(dir1.increment()(1) != dir2.increment()(1))
3748 : return false;
3749 : SpectralCoordinate sp1(coordsys.spectralCoordinate(2));
3750 : SpectralCoordinate sp2(imageCoord.spectralCoordinate(2));
3751 : if(sp1.increment()(0) != sp2.increment()(0))
3752 : return false;
3753 : */
3754 0 : return true;
3755 0 : }
3756 :
3757 0 : void Imager::setImageParam(Int& nx, Int& ny, Int& npol, Int& nchan){
3758 :
3759 0 : nx_p=nx;
3760 0 : ny_p=ny;
3761 0 : npol_p=npol;
3762 0 : nchan_p=nchan;
3763 :
3764 0 : }
3765 :
3766 53 : void Imager::makeVisSet(MeasurementSet& ms,
3767 : Bool compress, Bool mosaicOrder){
3768 :
3769 53 : if(rvi_p) {
3770 0 : delete rvi_p;
3771 0 : rvi_p=0;
3772 0 : wvi_p=0;
3773 : }
3774 :
3775 53 : Block<Int> sort(0);
3776 53 : if(mosaicOrder){
3777 0 : sort.resize(4);
3778 0 : sort[0] = MS::FIELD_ID;
3779 0 : sort[1] = MS::ARRAY_ID;
3780 0 : sort[2] = MS::DATA_DESC_ID;
3781 0 : sort[3] = MS::TIME;
3782 :
3783 : }
3784 : //else use default sort order
3785 : else{
3786 53 : sort.resize(4);
3787 53 : sort[0] = MS::ARRAY_ID;
3788 53 : sort[1] = MS::FIELD_ID;
3789 53 : sort[2] = MS::DATA_DESC_ID;
3790 53 : sort[3] = MS::TIME;
3791 : }
3792 53 : Matrix<Int> noselection;
3793 53 : Double timeInterval=0.0;
3794 : //if you want to use scratch col...make sure they are there
3795 53 : if(useModelCol_p){
3796 : //VisSet(ms,sort,noselection,useModelCol_p,timeInterval,compress);
3797 25 : VisSetUtil::addScrCols(ms, true, false, true, compress);
3798 : //delete keyword models to make sure data column is read
3799 25 : VisModelData::clearModel(ms);
3800 : }
3801 53 : if(imwgt_p.getType()=="none"){
3802 7 : imwgt_p=VisImagingWeight("natural");
3803 : }
3804 :
3805 53 : if(!ms.isWritable()){
3806 0 : rvi_p=new ROVisibilityIterator(ms, sort, timeInterval);
3807 : }
3808 : else{
3809 53 : wvi_p=new VisibilityIterator(ms, sort, timeInterval);
3810 53 : rvi_p=wvi_p;
3811 : }
3812 53 : rvi_p->useImagingWeight(imwgt_p);
3813 :
3814 : //////////////////////
3815 : //rvi_p->setRowBlocking(35);
3816 : ////////////////////
3817 53 : }
3818 : /*
3819 : void Imager::makeVisSet(MeasurementSet& ms,
3820 : Bool compress, Bool mosaicOrder){
3821 :
3822 :
3823 : Block<Int> sort(0);
3824 : if(mosaicOrder){
3825 : sort.resize(4);
3826 : sort[0] = MS::FIELD_ID;
3827 : sort[1] = MS::ARRAY_ID;
3828 : sort[2] = MS::DATA_DESC_ID;
3829 : sort[3] = MS::TIME;
3830 : }
3831 : //else use default sort order
3832 : else{
3833 :
3834 : sort.resize(4);
3835 : sort[0] = MS::ARRAY_ID;
3836 : sort[1] = MS::FIELD_ID;
3837 : sort[2] = MS::DATA_DESC_ID;
3838 : sort[3] = MS::TIME;
3839 : }
3840 : Matrix<Int> noselection;
3841 : Double timeInterval=0;
3842 :
3843 : VisSet vs(ms,sort,noselection,useModelCol_p,timeInterval,compress);
3844 :
3845 : }
3846 : */
3847 :
3848 44 : void Imager::writeHistory(LogIO& os){
3849 :
3850 88 : LogIO oslocal(LogOrigin("Imager", "writeHistory"));
3851 : try{
3852 :
3853 44 : os.postLocally();
3854 44 : if(!hist_p.null())
3855 44 : hist_p->addMessage(os);
3856 0 : }catch (AipsError x) {
3857 : oslocal << LogIO::SEVERE << "Caught exception: " << x.getMesg()
3858 0 : << LogIO::POST;
3859 0 : }
3860 44 : }
3861 :
3862 46 : void Imager::writeCommand(LogIO& os){
3863 :
3864 92 : LogIO oslocal(LogOrigin("Imager", "writeCommand"));
3865 : try{
3866 46 : os.postLocally();
3867 46 : if(!hist_p.null())
3868 46 : hist_p->cliCommand(os);
3869 0 : }catch (AipsError x) {
3870 : oslocal << LogIO::SEVERE << "Caught exception: " << x.getMesg()
3871 0 : << LogIO::POST;
3872 0 : }
3873 46 : }
3874 :
3875 0 : Bool Imager::makePBImage(ImageInterface<Float>& pbImage,
3876 : Bool useSymmetricBeam){
3877 :
3878 0 : LogIO os(LogOrigin("Imager", "makePBImage()", WHERE));
3879 0 : CoordinateSystem imageCoord=pbImage.coordinates();
3880 0 : Int spectralIndex=imageCoord.findCoordinate(Coordinate::SPECTRAL);
3881 : SpectralCoordinate
3882 0 : spectralCoord=imageCoord.spectralCoordinate(spectralIndex);
3883 0 : Vector<String> units(1); units = "Hz";
3884 0 : spectralCoord.setWorldAxisUnits(units);
3885 0 : Vector<Double> spectralWorld(1);
3886 0 : Vector<Double> spectralPixel(1);
3887 0 : spectralPixel(0) = 0;
3888 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
3889 0 : Double freq = spectralWorld(0);
3890 0 : Quantity qFreq( freq, "Hz" );
3891 0 : String telName=imageCoord.obsInfo().telescope();
3892 0 : if(telName=="UNKNOWN"){
3893 : os << LogIO::SEVERE << "Telescope encoded in image in not known "
3894 0 : << LogIO::POST;
3895 0 : return false;
3896 : }
3897 :
3898 :
3899 0 : PBMath myPB(telName, useSymmetricBeam, qFreq);
3900 0 : return makePBImage(myPB, pbImage);
3901 :
3902 0 : }
3903 :
3904 0 : Bool Imager::makePBImage(const CoordinateSystem& imageCoord,
3905 : const String& telescopeName,
3906 : const String& diskPBName,
3907 : Bool useSymmetricBeam, Double diam){
3908 :
3909 0 : LogIO os(LogOrigin("Imager", "makePBImage()", WHERE));
3910 0 : Int spectralIndex=imageCoord.findCoordinate(Coordinate::SPECTRAL);
3911 : SpectralCoordinate
3912 0 : spectralCoord=imageCoord.spectralCoordinate(spectralIndex);
3913 0 : Vector<String> units(1); units = "Hz";
3914 0 : spectralCoord.setWorldAxisUnits(units);
3915 0 : Vector<Double> spectralWorld(1);
3916 0 : Vector<Double> spectralPixel(1);
3917 0 : spectralPixel(0) = 0;
3918 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
3919 0 : Double freq = spectralWorld(0);
3920 0 : Quantity qFreq( freq, "Hz" );
3921 0 : String telName=telescopeName;
3922 0 : if(telName=="ALMA" && diam < 12.0)
3923 0 : telName="ACA";
3924 : //cerr << "Telescope Name is " << telName<< " qFreq " << qFreq << endl;
3925 : PBMath::CommonPB whichPB;
3926 0 : PBMath::enumerateCommonPB(telName, whichPB);
3927 0 : PBMath myPB;
3928 0 : if(whichPB!=PBMath::UNKNOWN && whichPB!=PBMath::NONE){
3929 :
3930 0 : myPB=PBMath(telName, useSymmetricBeam, qFreq);
3931 : }
3932 0 : else if(diam > 0.0){
3933 0 : myPB=PBMath(diam,useSymmetricBeam, qFreq);
3934 : }
3935 : else{
3936 : os << LogIO::WARN << "Telescope " << telName << " is not known\n "
3937 : << "Not making the PB image"
3938 0 : << LogIO::POST;
3939 0 : return false;
3940 : }
3941 0 : return makePBImage(imageCoord, myPB, diskPBName);
3942 :
3943 0 : }
3944 :
3945 0 : Bool Imager::makePBImage(const CoordinateSystem& imageCoord,
3946 : const Table& vpTable, const String& diskPBName){
3947 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
3948 0 : PBMath myPB(recCol(0));
3949 0 : return makePBImage(imageCoord, myPB, diskPBName);
3950 :
3951 0 : }
3952 :
3953 :
3954 0 : Bool Imager::makePBImage(const Table& vpTable, ImageInterface<Float>& pbImage){
3955 0 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
3956 0 : PBMath myPB(recCol(0));
3957 0 : return makePBImage(myPB, pbImage);
3958 :
3959 0 : }
3960 :
3961 0 : Bool Imager::makePBImage(const CoordinateSystem& /*imageCoord*/, PBMath& pbMath,
3962 : const String& diskPBName){
3963 :
3964 0 : make(diskPBName);
3965 0 : PagedImage<Float> pbImage(diskPBName);
3966 0 : return makePBImage(pbMath, pbImage);
3967 0 : }
3968 :
3969 :
3970 0 : Bool Imager::makePBImage(PBMath& pbMath, ImageInterface<Float>& pbImage){
3971 :
3972 0 : CoordinateSystem imageCoord=pbImage.coordinates();
3973 0 : pbImage.set(0.0);
3974 0 : MDirection wcenter;
3975 0 : Int directionIndex=imageCoord.findCoordinate(Coordinate::DIRECTION);
3976 : DirectionCoordinate
3977 0 : directionCoord=imageCoord.directionCoordinate(directionIndex);
3978 :
3979 0 : IPosition imShape=pbImage.shape();
3980 : //Vector<Double> pcenter(2);
3981 : // pcenter(0) = imShape(0)/2;
3982 : // pcenter(1) = imShape(1)/2;
3983 : //directionCoord.toWorld( wcenter, pcenter );
3984 0 : VisBuffer vb(*rvi_p);
3985 0 : Int fieldCounter=0;
3986 0 : Vector<Int> fieldsDone;
3987 :
3988 0 : for(rvi_p->originChunks(); rvi_p->moreChunks(); rvi_p->nextChunk()){
3989 0 : Bool fieldDone=false;
3990 0 : for (uInt k=0; k < fieldsDone.nelements(); ++k)
3991 : //hopefully there is not more than 10000 fields per ms
3992 0 : fieldDone=fieldDone || ((vb.msId()*10000)+vb.fieldId())==fieldsDone(k);
3993 0 : if(!fieldDone){
3994 0 : ++fieldCounter;
3995 0 : fieldsDone.resize(fieldCounter, true);
3996 0 : fieldsDone(fieldCounter-1)=vb.fieldId()+vb.msId()*10000;
3997 0 : wcenter=vb.msColumns().field().phaseDirMeas(vb.fieldId());
3998 0 : TempImage<Float> pbTemp(imShape, imageCoord);
3999 0 : TempImage<Complex> ctemp(imShape, imageCoord);
4000 0 : ctemp.set(1.0);
4001 0 : pbMath.applyPB(ctemp, ctemp, wcenter, Quantity(0.0, "deg"), BeamSquint::NONE);
4002 0 : StokesImageUtil::To(pbTemp, ctemp);
4003 0 : pbImage.copyData( (LatticeExpr<Float>)(pbImage+pbTemp) );
4004 0 : }
4005 : }
4006 0 : LatticeExprNode elmax= max( pbImage );
4007 0 : Float fmax = abs(elmax.getFloat());
4008 : //If there are multiple overlap of beam such that the peak is larger than 1 then normalize
4009 : //otherwise leave as is
4010 0 : if(fmax>1.0)
4011 0 : pbImage.copyData((LatticeExpr<Float>)(pbImage/fmax));
4012 :
4013 0 : Float cutoffval=minPB_p;
4014 0 : LatticeExpr<Bool> lemask(iif(pbImage < cutoffval,
4015 0 : false, true));
4016 0 : ImageRegion outreg=pbImage.makeMask("mask0", false, true);
4017 0 : LCRegion& outmask=outreg.asMask();
4018 0 : outmask.copyData(lemask);
4019 0 : pbImage.defineRegion("mask0", outreg,RegionHandler::Masks, true);
4020 0 : pbImage.setDefaultMask("mask0");
4021 :
4022 :
4023 :
4024 0 : return true;
4025 0 : }
4026 0 : void Imager::transferHistory(LoggerHolder& imagelog, MSHistoryColumns& msHis){
4027 0 : LogIO os(LogOrigin("imager", "transferHistory"));
4028 0 : LogSink& sink = imagelog.sink();
4029 0 : const ScalarColumn<Double> &time_col = msHis.time();
4030 0 : const ScalarColumn<String> &origin_col = msHis.origin();
4031 0 : const ArrayColumn<String> &cli_col = msHis.cliCommand();
4032 0 : const ScalarColumn<String> &message_col = msHis.message();
4033 0 : const ScalarColumn<String> &priority_col = msHis.priority();
4034 0 : std::map<String, LogMessage::Priority> prio={{"DEBUGGING", LogMessage::DEBUGGING},
4035 0 : {"DEBUG2", LogMessage::DEBUG2},{"DEBUG1", LogMessage::DEBUG1}, {"NORMAL5", LogMessage::NORMAL5}, {"NORMAL4", LogMessage::NORMAL4}, {"NORMAL3", LogMessage::NORMAL3}, {"NORMAL2", LogMessage::NORMAL2}, {"NORMAL1", LogMessage::NORMAL1}, {"NORMAL", LogMessage::NORMAL}, {"WARN", LogMessage::WARN}, {"SEVERE", LogMessage::SEVERE}};
4036 0 : if (msHis.nrow()>0) {
4037 : //ostringstream oos;
4038 0 : uInt nmessages = time_col.nrow();
4039 0 : for (uInt i=0; i < nmessages; i++) {
4040 : try{
4041 0 : ostringstream oos;
4042 0 : oos << cli_col(i);
4043 0 : LogMessage msg1(String("CLI_COMM: " + String(oos) + "; MESSAGE: " +message_col(i)), LogOrigin(origin_col(i)), prio.find(priority_col(i)) != prio.end() ? prio[priority_col(i)] :LogMessage::NORMAL);
4044 0 : msg1.messageTime( MVTime(Quantity(time_col(i), "s")).getTime());
4045 :
4046 : /* String tmp=frmtTime(time_col(i));
4047 : oos << tmp
4048 : << " HISTORY " << origin_col(i);
4049 : oos << " " << cli_col(i) << " ";
4050 : oos << message_col(i)
4051 : << endl;*/
4052 :
4053 0 : sink.postLocally(msg1);
4054 0 : }
4055 0 : catch(exception& y){
4056 0 : os << LogIO::DEBUG2 << "Skipping history-table row " << i << " while filling output logsink-header " << LogIO::POST;
4057 0 : }
4058 :
4059 : }
4060 : // String historyline(oos);
4061 : //sink.postLocally(msg.message(oos.str()));
4062 : }
4063 :
4064 0 : }
4065 0 : void Imager::setObsInfo(ObsInfo& obsinfo){
4066 :
4067 0 : latestObsInfo_p=obsinfo;
4068 0 : }
4069 :
4070 0 : ObsInfo& Imager::latestObsInfo(){
4071 0 : return latestObsInfo_p;
4072 : }
4073 :
4074 0 : Bool Imager::makeEmptyImage(CoordinateSystem& coords, String& name, Int fieldID){
4075 :
4076 0 : Int tilex=32;
4077 0 : if(imageTileVol_p >0){
4078 0 : tilex=static_cast<Int>(ceil(sqrt(imageTileVol_p/min(4, npol_p)/min(32, imageNchan_p))));
4079 0 : if(tilex >0){
4080 0 : if(tilex > min(nx_p, ny_p))
4081 0 : tilex=min(nx_p, ny_p);
4082 : else
4083 0 : tilex=nx_p/Int(nx_p/tilex);
4084 : }
4085 : //Not too small in x-y tile
4086 0 : if(tilex < 10)
4087 0 : tilex=10;
4088 : }
4089 0 : IPosition tileShape(4, min(tilex, nx_p), min(tilex, ny_p),
4090 0 : min(4, npol_p), min(32, imageNchan_p));
4091 0 : IPosition imageShape(4, nx_p, ny_p, npol_p, imageNchan_p);
4092 0 : PagedImage<Float> modelImage(TiledShape(imageShape, tileShape), coords, name);
4093 0 : modelImage.set(0.0);
4094 0 : modelImage.table().markForDelete();
4095 :
4096 : // Fill in miscellaneous information needed by FITS
4097 0 : MSColumns msc(*ms_p);
4098 0 : Record info;
4099 0 : String object=msc.field().name()(fieldID)
4100 : ; //defining object name
4101 0 : String objectName=msc.field().name()(fieldID);
4102 0 : ImageInfo iinfo=modelImage.imageInfo();
4103 0 : iinfo.setObjectName(objectName);
4104 0 : modelImage.setImageInfo(iinfo);
4105 0 : String telescop=msc.observation().telescopeName()(0);
4106 0 : info.define("INSTRUME", telescop);
4107 0 : info.define("distance", 0.0);
4108 0 : modelImage.setMiscInfo(info);
4109 0 : modelImage.table().tableInfo().setSubType("GENERIC");
4110 0 : modelImage.setUnits(Unit("Jy/beam"));
4111 0 : modelImage.table().unmarkForDelete();
4112 0 : modelImage.table().relinquishAutoLocks(true);
4113 0 : modelImage.table().unlock();
4114 :
4115 0 : return true;
4116 :
4117 0 : }
4118 :
4119 0 : String Imager::frmtTime(const Double time) {
4120 0 : MVTime mvtime(Quantity(time, "s"));
4121 0 : Time t=mvtime.getTime();
4122 0 : ostringstream os;
4123 0 : os << t;
4124 0 : return os.str();
4125 0 : }
4126 :
4127 0 : Bool Imager::getRestFreq(Vector<Double>& restFreq, const Int& spw){
4128 : // MS Doppler tracking utility
4129 0 : MSDopplerUtil msdoppler(*ms_p);
4130 0 : restFreq.resize();
4131 0 : if(restFreq_p.getValue() > 0){// User defined restfrequency
4132 0 : restFreq.resize(1);
4133 0 : restFreq[0]=restFreq_p.getValue("Hz");
4134 : }
4135 : else{
4136 : // Look up first rest frequency found (for now)
4137 :
4138 0 : Int fieldid = (datafieldids_p.nelements()>0 ? datafieldids_p(0) :
4139 0 : fieldid_p);
4140 : try{
4141 0 : msdoppler.dopplerInfo(restFreq ,spw,fieldid);
4142 : }
4143 0 : catch(...){
4144 0 : restFreq.resize();
4145 0 : }
4146 : }
4147 0 : if(restFreq.nelements() >0)
4148 0 : return true;
4149 0 : return false;
4150 0 : }
4151 :
4152 :
4153 37 : void Imager::setSkyEquation(){
4154 : /* if(sm_p->nmodels() >0){
4155 : Long npix=0;
4156 : for (Int model=0; model < sm_p->nmodels(); ++model){
4157 : Long pixmod=sm_p->image(model).product();
4158 : npix=max(pixmod, npix);
4159 : }
4160 : Long pixInMem=(HostInfo::memoryTotal()/8)*1024;
4161 : if(npix > (pixInMem/2)){
4162 : se_p = new CubeSkyEquation(*sm_p, *vs_p, *ft_p, *cft_p, !useModelCol_p);
4163 : return;
4164 : }
4165 : //else lets make the base SkyEquation for now
4166 : }
4167 :
4168 : se_p = new SkyEquation(*sm_p, *vs_p, *ft_p, *cft_p, !useModelCol_p);
4169 :
4170 : */
4171 : // if (ft_p->name() == "PBWProjectFT")
4172 : // {
4173 : // logSink_p.clearLocally();
4174 : // LogIO os(LogOrigin("imager", "setSkyEquation()"), logSink_p);
4175 : // os << "Creating SkyEquation for PBWProjectFT" << LogIO::POST;
4176 : // se_p = new SkyEquation(*sm_p, *vs_p, *ft_p, *cft_p, !useModelCol_p);
4177 : // }
4178 : // else
4179 37 : se_p = new CubeSkyEquation(*sm_p, *rvi_p, *ft_p, *cft_p, !useModelCol_p);
4180 37 : return;
4181 : }
4182 :
4183 0 : void Imager::savePSF(const Vector<String>& psf){
4184 :
4185 0 : if( (psf.nelements() > 0) && (Int(psf.nelements()) <= sm_p->numberOfModels())){
4186 :
4187 0 : for (Int thismodel=0;thismodel<Int(psf.nelements());++thismodel) {
4188 0 : if(removeTable(psf(thismodel))) {
4189 0 : Int whichmodel=thismodel;
4190 0 : if(facets_p >1 && thismodel > 0)
4191 0 : whichmodel=facets_p*facets_p-1+thismodel;
4192 0 : IPosition shape=images_p[thismodel]->shape();
4193 : PagedImage<Float> psfimage(shape,
4194 0 : images_p[thismodel]->coordinates(),
4195 0 : psf(thismodel));
4196 0 : if(freqFrameValid_p){
4197 0 : CoordinateSystem cs=psfimage.coordinates();
4198 0 : String errorMsg;
4199 0 : if (cs.setSpectralConversion (errorMsg, MFrequency::showType(freqFrame_p))) {
4200 0 : psfimage.setCoordinateInfo(cs);
4201 : }
4202 0 : }
4203 0 : psfimage.set(0.0);
4204 0 : if((shape[0]*shape[1]) > ((sm_p->PSF(whichmodel)).shape()[0]*(sm_p->PSF(whichmodel)).shape()[1])){
4205 0 : IPosition blc(4, 0, 0, 0, 0);
4206 0 : IPosition trc=shape-1;
4207 0 : blc[0]=(shape[0]-(sm_p->PSF(whichmodel)).shape()[0])/2;
4208 0 : trc[0]=(sm_p->PSF(whichmodel)).shape()[0]+blc[0]-1;
4209 0 : blc[1]=(shape[1]-(sm_p->PSF(whichmodel)).shape()[1])/2;
4210 0 : trc[1]=(sm_p->PSF(whichmodel)).shape()[1]+blc[1]-1;
4211 0 : Slicer sl(blc, trc, Slicer::endIsLast);
4212 0 : SubImage<Float> sub(psfimage, sl, true);
4213 0 : sub.copyData(sm_p->PSF(whichmodel));
4214 0 : }
4215 : else{
4216 0 : psfimage.copyData(sm_p->PSF(whichmodel));
4217 : }
4218 : //sm_p->PSF(whichmodel).clearCache();
4219 0 : }
4220 : }
4221 :
4222 :
4223 :
4224 : }
4225 :
4226 0 : }
4227 :
4228 0 : void Imager::setMosaicFTMachine(Bool useDoublePrec){
4229 0 : LogIO os(LogOrigin("Imager", "setmosaicftmachine", WHERE));
4230 0 : MSColumns msc(*ms_p);
4231 0 : String telescop=msc.observation().telescopeName()(0);
4232 : PBMath::CommonPB kpb;
4233 0 : PBMath::enumerateCommonPB(telescop, kpb);
4234 0 : if(!((kpb == PBMath::UNKNOWN) ||
4235 0 : (kpb==PBMath::OVRO) || (kpb==PBMath::ALMA) || (kpb==PBMath::ACA) || (kpb==PBMath::EVLA) || !doDefaultVP_p)){
4236 :
4237 0 : if(!gvp_p) {
4238 0 : MSColumns msc(*ms_p);
4239 0 : if (doDefaultVP_p) {
4240 : os << LogIO::NORMAL // Loglevel INFO
4241 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
4242 0 : gvp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p,
4243 0 : skyPosThreshold_p);
4244 : } /*else {
4245 : os << LogIO::NORMAL // Loglevel INFO
4246 : << "Using VP as defined in " << vpTableStr_p << LogIO::POST;
4247 : Table vpTable( vpTableStr_p );
4248 : gvp_p=new VPSkyJones(msc, vpTable, parAngleInc_p, squintType_p,
4249 : skyPosThreshold_p);
4250 : }
4251 : */
4252 0 : }
4253 0 : gvp_p->setThreshold(minPB_p);
4254 : }
4255 :
4256 0 : ft_p = new MosaicFT(gvp_p, mLocation_p, stokes_p, cache_p/2, tile_p, true,
4257 0 : useDoublePrec);
4258 :
4259 0 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
4260 0 : || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || (!doDefaultVP_p)){
4261 : os << LogIO::NORMAL // Loglevel INFO
4262 : << "Using antenna primary beams for determining beams for gridding"
4263 0 : << LogIO::POST;
4264 : //Use 1D-Airy
4265 0 : PBMathInterface::PBClass pbtype= (kpb==PBMath::EVLA) ? PBMathInterface::COMMONPB : PBMathInterface::AIRY;
4266 : //Temporary special case for ALMA to use 2D images
4267 : //Temporary till it is made fully that with automatic image selections
4268 0 : if((kpb==PBMath::ACA) || (kpb==PBMath::ALMA)){
4269 0 : String useimage="false";
4270 0 : Aipsrc::find(useimage, "alma.vp.image", "false");
4271 0 : useimage.downcase();
4272 0 : if(useimage.contains("true")){
4273 0 : pbtype=PBMathInterface::IMAGE;
4274 :
4275 : }
4276 0 : }
4277 0 : if(!doDefaultVP_p){
4278 0 : pbtype=PBMathInterface::IMAGE;
4279 : }
4280 : else{
4281 0 : vpTableStr_p="";
4282 : }
4283 0 : CountedPtr<SimplePBConvFunc> mospb=new HetArrayConvFunc(pbtype, vpTableStr_p);
4284 :
4285 0 : static_cast<MosaicFT &>(*ft_p).setConvFunc(mospb);
4286 0 : }
4287 :
4288 0 : }
4289 0 : ATerm* Imager::createTelescopeATerm(MeasurementSet& ms, const Bool& isATermOn)
4290 : {
4291 0 : LogIO log_l(LogOrigin("Imager", "createTelescopeATerm"));
4292 :
4293 0 : if (!isATermOn) return new NoOpATerm();
4294 :
4295 0 : MSObservationColumns msoc(ms.observation());
4296 0 : String ObsName=msoc.telescopeName()(0);
4297 0 : if ((ObsName == "EVLA") || (ObsName == "VLA"))
4298 0 : return new EVLAAperture();
4299 : else
4300 : {
4301 0 : log_l << "Telescope name ('"+
4302 0 : ObsName+"') in the MS not recognized to create the telescope specific ATerm"
4303 0 : << LogIO::WARN;
4304 : }
4305 :
4306 0 : return NULL;
4307 0 : }
4308 :
4309 : //use SubMS::calcChanFreqs to calculate spectral gridding
4310 : //call from imagecoodinates2
4311 0 : Bool Imager::calcImFreqs(Vector<Double>& imgridfreqs,
4312 : Vector<Double>& imfreqres,
4313 : const MFrequency::Types& oldRefFrame,
4314 : const MEpoch& obsEpoch,
4315 : const MPosition& obsPosition,
4316 : const Double& restFreq
4317 : )
4318 : {
4319 :
4320 0 : logSink_p.clearLocally();
4321 0 : LogIO os(LogOrigin("imager", "setGridFreqs()"), logSink_p);
4322 :
4323 0 : MSColumns msc(*ms_p);
4324 0 : Vector<Double> oldChanFreqs;
4325 0 : Vector<Double> oldFreqResolution;
4326 0 : String veltype;
4327 0 : String mode;
4328 0 : String restfreq;
4329 0 : String start;
4330 0 : String width;
4331 0 : String outframe;
4332 0 : Bool reversevec(false);
4333 0 : Bool descendfreq(false);
4334 :
4335 0 : if (imageMode_p.contains("RADIO")) {
4336 0 : veltype="radio";
4337 0 : mode="velocity";
4338 0 : start=dQuantitytoString(mImageStart_p.get("m/s"));
4339 0 : width=dQuantitytoString(mImageStep_p.get("m/s"));
4340 0 : if (!width.contains(casacore::Regex("^-"))) {
4341 : //positive vel. width (descending frequencies)
4342 : //reversevec=true;
4343 0 : descendfreq=true;
4344 : }
4345 : }
4346 0 : else if (imageMode_p.contains("OPTICAL")) {
4347 0 : veltype="optical";
4348 0 : mode="velocity";
4349 0 : start=dQuantitytoString(mImageStart_p.get("m/s"));
4350 0 : width=dQuantitytoString(mImageStep_p.get("m/s"));
4351 : //cerr<<"optical vel width USED="<<width<<endl;
4352 0 : if (!width.contains(casacore::Regex("^-"))) {
4353 : //positive vel. width (descending frequencies)
4354 : //reversevec=true;
4355 0 : descendfreq=true;
4356 : }
4357 : }
4358 0 : else if (imageMode_p.contains("FREQ")) {
4359 0 : veltype="radio";
4360 0 : mode="frequency";
4361 0 : start=dQuantitytoString(mfImageStart_p.get("Hz"));
4362 0 : width=dQuantitytoString(mfImageStep_p.get("Hz"));
4363 0 : if (width.contains(casacore::Regex("^-"))) {
4364 : //reversevec=true;
4365 0 : descendfreq=true;
4366 : }
4367 : }
4368 0 : else if (imageMode_p.contains("CHANNEL")) {
4369 0 : veltype="radio";
4370 0 : mode="channel";
4371 0 : start=String::toString(imageStart_p);
4372 0 : width=String::toString(imageStep_p);
4373 0 : if (width.contains(casacore::Regex("^-"))) {
4374 : // means here going to lower chan index
4375 0 : descendfreq=true;
4376 : }
4377 : }
4378 0 : else if (imageMode_p.contains("MFS")) {
4379 0 : veltype="radio";
4380 0 : mode="mfs";
4381 0 : start=String::toString(imageStart_p);
4382 : }
4383 :
4384 0 : restfreq = dQuantitytoString(Quantity(restFreq,"Hz"));
4385 : //MFrequency::getType(freqFrame_p, outframe);
4386 0 : if (freqFrame_p!=oldRefFrame) {
4387 0 : outframe=MFrequency::showType(freqFrame_p);
4388 : }
4389 : else {
4390 0 : outframe="";
4391 : }
4392 :
4393 :
4394 : try {
4395 : /***
4396 : if(spectralwindowids_p.nelements()==1){
4397 : if(spectralwindowids_p[0]<0){
4398 : spectralwindowids_p.resize();
4399 : if(dataspectralwindowids_p.nelements()==0){
4400 : Int nspwinms=ms_p->spectralWindow().nrow();
4401 : dataspectralwindowids_p.resize(nspwinms);
4402 : indgen(dataspectralwindowids_p);
4403 : }
4404 : spectralwindowids_p=dataspectralwindowids_p;
4405 : }
4406 : }
4407 : ***/
4408 0 : Vector<Int> spwlist;
4409 0 : if (mode=="frequency" || mode=="velocity") {
4410 0 : spwlist = dataspectralwindowids_p;
4411 : }
4412 : else {
4413 0 : spwlist = spectralwindowids_p;
4414 : }
4415 0 : if(spwlist.nelements()==1) {
4416 0 : oldChanFreqs=msc.spectralWindow().chanFreq()(spwlist[0]);
4417 0 : oldFreqResolution=msc.spectralWindow().chanWidth()(spwlist[0]);
4418 : }
4419 : else {
4420 0 : SubMS thems(*ms_p);
4421 0 : if(!thems.combineSpws(spwlist,true,oldChanFreqs,oldFreqResolution)){
4422 0 : os << LogIO::SEVERE << "Error combining SpWs" << LogIO::POST;
4423 : }
4424 0 : }
4425 0 : Bool isDescendingData=false;
4426 : // Some descending order data has positive channel widths...so check chan freqs
4427 : // first...
4428 : //if (oldFreqResolution(0) < 0) isDescendingData=true;
4429 0 : if (oldChanFreqs.nelements()>1) {
4430 0 : if ((oldChanFreqs[1] - oldChanFreqs[0])<0) isDescendingData=true;
4431 : }
4432 0 : else if (oldFreqResolution(0) < 0) {
4433 0 : isDescendingData=true;
4434 : }
4435 :
4436 : // need theOldRefFrame,theObsTime,mObsPos,mode,nchan,start,width,restfreq,
4437 : // outframe,veltype
4438 : //
4439 :
4440 0 : Bool rst=SubMS::calcChanFreqs(os,
4441 : imgridfreqs,
4442 : imfreqres,
4443 : oldChanFreqs,
4444 : oldFreqResolution,
4445 0 : phaseCenter_p,
4446 : oldRefFrame,
4447 : obsEpoch,
4448 : obsPosition,
4449 : mode,
4450 : imageNchan_p,
4451 : start,
4452 : width,
4453 : restfreq,
4454 : outframe,
4455 : veltype
4456 : );
4457 :
4458 0 : if (!rst) {
4459 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
4460 0 : << LogIO::EXCEPTION;
4461 0 : return false;
4462 : }
4463 :
4464 : //cout<<"=================calcChanFreqs arguments ======================"<<endl;
4465 : //cout.precision(10);
4466 : //cout<<"imgridfreqs(0)="<<imgridfreqs(0)<<endl;
4467 : //cout<<"imgridfreqs("<<imgridfreqs.nelements()-1<<")="<<imgridfreqs(imgridfreqs.nelements()-1)<<endl;
4468 : //cout<<"oldChanFreqs(0)="<<oldChanFreqs(0)<<endl;
4469 : //cout<<"oldChanFreqs("<<oldChanFreqs.nelements()-1<<")="<<oldChanFreqs(oldChanFreqs.nelements()-1)<<endl;
4470 : //cout<<"phaseCenter_p="<<phaseCenter_p<<endl;
4471 : //cout<<"oldRefFrame="<<oldRefFrame<<endl;
4472 : //cout<<"outframe="<<outframe<<endl;
4473 : //cout<<"obsEpoch="<<obsEpoch<<endl;
4474 : //cout<<"obsPosition="<<obsPosition<<endl;
4475 : //cout<<"start="<<start<<" width="<<width<<endl;
4476 : //cout<<"restfreq="<<restfreq<<endl;
4477 : //cout<<"veltype="<<veltype<<endl;
4478 : //cout<<"=================calcChanFreqs arguments end==================="<<endl;
4479 0 : Bool isDescendingNewData=false;
4480 0 : if (imgridfreqs(0)-imgridfreqs(1)>0) isDescendingNewData=true;
4481 : //reverse frequency vector?
4482 : //evaluate reversing condition differ for chan mode from other modes
4483 0 : if(mode.contains("channel")) {
4484 0 : if ((descendfreq && !isDescendingNewData && !isDescendingData) |
4485 0 : (descendfreq && isDescendingNewData && isDescendingData) |
4486 0 : (!descendfreq && !isDescendingNewData && isDescendingData)) {
4487 0 : reversevec = true;
4488 : }
4489 : }
4490 : else {
4491 0 : if ((descendfreq && !isDescendingNewData) |
4492 0 : (!descendfreq && isDescendingNewData)) {
4493 0 : reversevec = true;
4494 : }
4495 : }
4496 0 : if (reversevec) {
4497 : //if(reversevec && isAscendingData ) {
4498 : //Int ndata=imgridfreqs.nelements();
4499 : //tempimgridfreqs.resize(ndata);
4500 : /**
4501 : for (Int i=0;i<ndata;i++) {
4502 : tempimgridfreqs[i] = imgridfreqs[ndata - 1 - i];
4503 : }
4504 : for (Int i=0;i<ndata;i++) {
4505 : std::swap(imgridfreqs[ndata-1-i],tempimgridfreqs[i]);
4506 : }
4507 : **/
4508 0 : std::vector<double> stlimgridfreqs;
4509 0 : imgridfreqs.tovector(stlimgridfreqs);
4510 0 : std::reverse(stlimgridfreqs.begin(),stlimgridfreqs.end());
4511 0 : imgridfreqs=Vector<Double>(stlimgridfreqs);
4512 0 : }
4513 : //cerr<<"Final imgridfreqs(0)="<<imgridfreqs(0)<<endl;
4514 :
4515 0 : } catch (AipsError x) {
4516 0 : this->unlock();
4517 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4518 0 : << LogIO::EXCEPTION;
4519 0 : return false;
4520 0 : }
4521 0 : return true;
4522 0 : }//end of calcImFreqs
4523 :
4524 : // convert a double precision quanity to a String
4525 0 : String Imager::dQuantitytoString(const Quantity& dq) {
4526 0 : std::ostringstream ss;
4527 0 : ss.precision(std::numeric_limits<double>::digits10);
4528 0 : ss << dq;
4529 0 : return ss.str();
4530 0 : }
4531 :
4532 : #define INITIALIZE_DIRECTION_VECTOR(name) \
4533 : do { \
4534 : if (name.nelements() < 2) { \
4535 : name.resize(2); \
4536 : name = 0.0; \
4537 : } \
4538 : } while (false)
4539 :
4540 0 : Bool Imager::mapExtent(const String &referenceFrame, const String &movingSource,
4541 : const String &pointingColumn, Vector<Double> ¢er, Vector<Double> &blc,
4542 : Vector<Double> &trc, Vector<Double> &extent) {
4543 0 : return getMapExtent(*mssel_p, referenceFrame, movingSource, pointingColumn,
4544 0 : center, blc, trc, extent);
4545 : }
4546 :
4547 0 : Bool Imager::getMapExtent(const MeasurementSet &ms, const String &referenceFrame,
4548 : const String &movingSource, const String &pointingColumn,
4549 : Vector<Double> ¢er, Vector<Double> &blc, Vector<Double> &trc, Vector<Double> &extent) {
4550 0 : INITIALIZE_DIRECTION_VECTOR(center);
4551 0 : INITIALIZE_DIRECTION_VECTOR(blc);
4552 0 : INITIALIZE_DIRECTION_VECTOR(trc);
4553 0 : INITIALIZE_DIRECTION_VECTOR(extent);
4554 :
4555 : try {
4556 : //cout << "ms.nrow() = " << ms.nrow() << endl;
4557 0 : PointingDirectionCalculator calc(ms);
4558 : //cout << "calc instantiated" << endl;
4559 :
4560 0 : calc.setDirectionColumn(pointingColumn);
4561 : //cout << "set pointing direction column to " << pointingColumn << endl;
4562 0 : calc.setFrame(referenceFrame);
4563 : //cout << "set reference frame to " << referenceFrame << endl;
4564 0 : MDirection::Types refType = MDirection::J2000; // any non planet value
4565 0 : Bool status = false;
4566 : //cout << "movingSource = " << movingSource << endl;
4567 : //cout << "getType" << endl;
4568 0 : status = MDirection::getType(refType, movingSource);
4569 : //cout << "refType string = " << MDirection::showType(refType) << endl;
4570 : //cout << "status = " << status << endl;
4571 0 : Bool doMovingSourceCorrection = (status == true &&
4572 0 : MDirection::N_Types < refType &&
4573 0 : refType < MDirection::N_Planets);
4574 0 : Bool isOffsetColumn = (pointingColumn.contains("OFFSET")
4575 0 : || pointingColumn.contains("Offset")
4576 0 : || pointingColumn.contains("offset"));
4577 0 : if (doMovingSourceCorrection) {
4578 : //cout << "need moving source correction" << endl;
4579 0 : calc.setMovingSource(movingSource);
4580 : //cout << "set moving source name " << movingSource << endl;
4581 : }
4582 : //cout << "set direction matrix shape to COLUMN_MAJOR" << endl;
4583 0 : calc.setDirectionListMatrixShape(PointingDirectionCalculator::COLUMN_MAJOR);
4584 :
4585 : //cout << "start getDirection" << endl;
4586 0 : Matrix<Double> directionList = calc.getDirection();
4587 : //cout << "end getDirection shape " << directionList.shape() << endl;
4588 0 : Vector<Double> longitude = directionList.column(0);
4589 0 : Vector<Double> latitude = directionList.column(1);
4590 : //cout << "longitude.nelements() = " << longitude.nelements() << endl;
4591 : //cout << "latitude.nelements() = " << latitude.nelements() << endl;
4592 :
4593 : // Diagnose if longitude values are divided by periodic boundary surface
4594 : // (+-pi or 0, 2pi)
4595 : // In this case, mean of longitude should be around 0 (+-pi) or pi (0,2pi)
4596 : // and stddev of longitude array be around pi.
4597 : //cout << "diagnose longitude distribution" << endl;
4598 0 : Double longitudeMean = mean(longitude);
4599 0 : Double longitudeStddev = stddev(longitude);
4600 : //cout << "mean " << longitudeMean << " stddev " << longitudeStddev << endl;
4601 0 : if (longitudeStddev > 2.0 / 3.0 * C::pi) {
4602 : // likely to be the case
4603 : //cout << "likely to be the case" << endl;
4604 0 : if (abs(longitudeMean) < 0.5 * C::pi) {
4605 : // periodic boundary surface would be +-pi
4606 : //cout << "periodic boundary surface would be +-pi" << endl;
4607 0 : for (size_t i = 0; i < longitude.nelements(); ++i) {
4608 0 : if (longitude[i] < 0.0) {
4609 0 : longitude[i] += C::_2pi;
4610 : }
4611 : }
4612 : }
4613 0 : else if (abs(longitudeMean - C::pi) < 0.5 * C::pi ) {
4614 : // periodic boundary surface would be 0,2pi
4615 : //cout << "periodic boundary surface would be 0,2pi" << endl;
4616 0 : for (size_t i = 0; i < longitude.nelements(); ++i) {
4617 0 : if (longitude[i] < C::pi) {
4618 0 : longitude[i] += C::_2pi;
4619 : }
4620 : }
4621 : }
4622 : }
4623 :
4624 0 : blc[0] = min(longitude);
4625 0 : trc[0] = max(longitude);
4626 0 : blc[1] = min(latitude);
4627 0 : trc[1] = max(latitude);
4628 0 : extent = trc - blc;
4629 : //cout << "result: blc " << blc << " trc " << trc << endl;
4630 : //cout << "result: extent " << extent << endl;
4631 :
4632 : // declination correction factor
4633 0 : Double declinationCorrection = 0.0;
4634 : // center
4635 0 : if (isOffsetColumn) {
4636 0 : center = 0.0;
4637 0 : declinationCorrection = cos((blc[1] + trc[1]) / 2.0);
4638 : }
4639 0 : else if (doMovingSourceCorrection) {
4640 : // shift center to moving source position at the time
4641 : // that will be used for imaging
4642 0 : VisBuffer vb(*rvi_p);
4643 0 : MEpoch referenceEpoch = vb.msColumns().timeMeas()(vb.rowIds()[0]);
4644 0 : MPosition referencePosition = vb.msColumns().antenna().positionMeas()(vb.antenna1()[0]);
4645 0 : MDirection::Types const &outDirectionType = calc.getDirectionType();
4646 0 : MDirection const &movingSourceDir = calc.getMovingSourceDirection();
4647 0 : MeasFrame referenceMeasFrame(referenceEpoch, referencePosition);
4648 0 : MDirection::Ref azelRef(MDirection::AZEL, referenceMeasFrame);
4649 0 : MDirection azelDir = MDirection::Convert(movingSourceDir, azelRef)();
4650 0 : MDirection::Ref outRef(outDirectionType, referenceMeasFrame);
4651 0 : MDirection sourceDirection = MDirection::Convert(azelDir, outRef)();
4652 0 : center = sourceDirection.getAngle("rad").getValue();
4653 0 : declinationCorrection = cos(center[1]);
4654 0 : }
4655 : else {
4656 0 : center = (blc + trc) / 2.0;
4657 0 : declinationCorrection = cos(center[1]);
4658 : }
4659 0 : assert(declinationCorrection != 0.0);
4660 : //cout << "declinationCorrection = " << declinationCorrection << endl;
4661 : //cout << "result: center " << center << endl;
4662 :
4663 : // apply declinationCorrection to extent[0]
4664 0 : extent[0] *= declinationCorrection;
4665 0 : }
4666 0 : catch (AipsError &e) {
4667 0 : LogIO os(LogOrigin("Imager", "getMapExtent", WHERE));
4668 0 : os << LogIO::SEVERE << "Failed due to the rror \"" << e.getMesg() << "\"" << LogIO::POST;
4669 0 : return false;
4670 0 : }
4671 0 : catch (...) {
4672 0 : LogIO os(LogOrigin("Imager", "getMapExtent", WHERE));
4673 0 : os << LogIO::SEVERE << "Failed due to unknown error" << LogIO::POST;
4674 0 : throw;
4675 : return false;
4676 0 : }
4677 :
4678 0 : return true;
4679 : }
4680 :
4681 0 : Bool Imager::pointingSampling(const String &referenceFrame,
4682 : const String &movingSource,
4683 : const String &pointingColumn,
4684 : const String &antenna,
4685 : Quantum<Vector<Double>> &sampling,
4686 : Quantity &positionAngle) {
4687 : try {
4688 0 : SingleDishBeamUtil sdBeamU(*mssel_p, referenceFrame, movingSource,
4689 0 : pointingColumn, antenna);
4690 0 : sdBeamU.getPointingSamplingRaster(sampling, positionAngle);
4691 0 : }
4692 0 : catch (AipsError &e) {
4693 0 : LogIO os(LogOrigin("Imager", "pointingSampling", WHERE));
4694 0 : os << LogIO::SEVERE << "Failed due to the rror \"" << e.getMesg() << "\"" << LogIO::POST;
4695 0 : return false;
4696 0 : }
4697 0 : catch (...) {
4698 0 : LogIO os(LogOrigin("Imager", "pointingSampling", WHERE));
4699 0 : os << LogIO::SEVERE << "Failed due to unknown error" << LogIO::POST;
4700 0 : throw;
4701 : return false;
4702 0 : }
4703 0 : return true;
4704 : }
4705 :
4706 : } //# NAMESPACE CASA - END
4707 :
4708 :
4709 :
4710 :
4711 :
4712 :
4713 : // else if ((ftmachine_p == "wbawp") || (ftmachine_p == "wbmosaic")){
4714 :
4715 : // if (wprojPlanes_p<=1)
4716 : // {
4717 : // os << LogIO::NORMAL
4718 : // << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
4719 : // << LogIO::POST;
4720 : // os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
4721 : // }
4722 : // if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
4723 : // {
4724 : // os << LogIO::WARN
4725 : // << "No. of w-planes set too low for W projection - recommend at least 128"
4726 : // << LogIO::POST;
4727 : // os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
4728 : // }
4729 :
4730 : // // if(!gvp_p)
4731 : // // {
4732 : // // os << LogIO::NORMAL // Loglevel INFO
4733 : // // << "Using defaults for primary beams used in gridding" << LogIO::POST;
4734 : // // gvp_p = new VPSkyJones(*ms_p, true, parAngleInc_p, squintType_p,
4735 : // // skyPosThreshold_p);
4736 : // // }
4737 : // useDoublePrecGrid=false;
4738 : // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(*ms_p,true);
4739 : // CountedPtr<PSTerm> psTerm = new PSTerm();
4740 : // CountedPtr<WTerm> wTerm = new WTerm();
4741 : // // psTerm->setOpCode(CFTerms::NOOP);
4742 : // CountedPtr<ConvolutionFunction> awConvFunc;
4743 : // if (ftmachine_p=="wbawp") awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm);
4744 : // else awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm);
4745 :
4746 : // CountedPtr<VisibilityResamplerBase> visResampler = new AWVisResampler();
4747 :
4748 : // // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
4749 :
4750 : // // CountedPtr<VisibilityResamplerBase> mthVisResampler =
4751 : // // new MultiThreadedVisibilityResampler(useDoublePrecGrid, visResampler);
4752 : // CountedPtr<CFCache> cfcache = new CFCache();
4753 : // cfcache->setCacheDir(cfCacheDirName_p.data());
4754 : // cerr << "cfcache->initCache2()" << endl;
4755 : // cfcache->initCache2();
4756 :
4757 : // ft_p = new AWProjectWBFT(wprojPlanes_p, cache_p/2,
4758 : // cfcache, awConvFunc,
4759 : // // mthVisResampler,
4760 : // visResampler,
4761 : // /*true */doPointing, doPBCorr,
4762 : // tile_p, paStep_p, pbLimit_p, true);
4763 :
4764 : // ((AWProjectWBFT *)ft_p)->setObservatoryLocation(mLocation_p);
4765 : // //
4766 : // // Explicit type casting of ft_p does not look good. It does not
4767 : // // pick up the setPAIncrement() method of PBWProjectFT without
4768 : // // this
4769 : // //
4770 : // // os << LogIO::NORMAL << "Setting PA increment to " << parAngleInc_p.getValue("deg") << " deg" << endl;
4771 : // ((AWProjectFT *)ft_p)->setPAIncrement(parAngleInc_p);
4772 :
4773 : // if (doPointing)
4774 : // {
4775 : // try
4776 : // {
4777 : // // Warn users we are have temporarily disabled pointing cal
4778 : // // throw(AipsError("Pointing calibration temporarily disabled (gmoellen 06Nov20)."));
4779 : // // TBD: Bring this up-to-date with new VisCal mechanisms
4780 : // VisSet elVS(*rvi_p);
4781 : // epJ = new EPJones(elVS, *ms_p);
4782 : // RecordDesc applyRecDesc;
4783 : // applyRecDesc.addField("table", TpString);
4784 : // applyRecDesc.addField("interp",TpString);
4785 : // Record applyRec(applyRecDesc);
4786 : // applyRec.define("table",epJTableName_p);
4787 : // applyRec.define("interp", "nearest");
4788 : // epJ->setApply(applyRec);
4789 : // ((AWProjectFT *)ft_p)->setEPJones(epJ);
4790 : // }
4791 : // catch(AipsError& x)
4792 : // {
4793 : // //
4794 : // // Add some more useful info. to the message and translate
4795 : // // the generic AipsError exception object to a more specific
4796 : // // SynthesisError object.
4797 : // //
4798 : // String mesg = x.getMesg();
4799 : // mesg += ". Error in loading pointing offset table.";
4800 : // SynthesisError err(mesg);
4801 : // throw(err);
4802 : // }
4803 : // }
4804 : // AlwaysAssert(ft_p, AipsError);
4805 : // cft_p = new SimpleComponentFTMachine();
4806 : // AlwaysAssert(cft_p, AipsError);
4807 :
4808 : // }
4809 : // else if (ftmachine_p == "awp")
4810 : // {
4811 : // if (wprojPlanes_p<=1)
4812 : // {
4813 : // os << LogIO::NORMAL
4814 : // << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
4815 : // << LogIO::POST;
4816 : // os << LogIO::NORMAL << "Performing A-Projection" << LogIO::POST; // Loglevel PROGRESS
4817 : // }
4818 : // if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
4819 : // {
4820 : // os << LogIO::WARN
4821 : // << "No. of w-planes set too low for W projection - recommend at least 128"
4822 : // << LogIO::POST;
4823 : // os << LogIO::NORMAL << "Performing AW-Projection"
4824 : // << LogIO::POST; // Loglevel PROGRESS
4825 : // }
4826 : // // if(!gvp_p)
4827 : // // {
4828 : // // os << LogIO::NORMAL // Loglevel INFO
4829 : // // << "Using defaults for primary beams used in gridding" << LogIO::POST;
4830 : // // gvp_p = new VPSkyJones(*ms_p, true, parAngleInc_p, squintType_p,
4831 : // // skyPosThreshold_p);
4832 : // // }
4833 : // // CountedPtr<ATerm> evlaAperture = new EVLAAperture();
4834 : // useDoublePrecGrid=false;
4835 : // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(*ms_p,true);
4836 : // CountedPtr<PSTerm> psTerm = new PSTerm();
4837 : // CountedPtr<WTerm> wTerm = new WTerm();
4838 : // psTerm->setOpCode(CFTerms::NOOP);
4839 : // CountedPtr<ConvolutionFunction> awConvFunc=new AWConvFunc(apertureFunction,psTerm,wTerm);
4840 : // CountedPtr<VisibilityResamplerBase> visResampler = new AWVisResampler();
4841 : // // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
4842 : // // CountedPtr<VisibilityResamplerBase> mthVisResampler = new MultiThreadedVisibilityResampler(useDoublePrecGrid,
4843 : // // visResampler);
4844 : // CountedPtr<CFCache> cfcache=new CFCache();
4845 : // cfcache->setCacheDir(cfCacheDirName_p.data());
4846 : // cfcache->initCache2();
4847 : // ft_p = new AWProjectFT(wprojPlanes_p, cache_p/2,
4848 : // cfcache, awConvFunc,
4849 : // // mthVisResampler,
4850 : // visResampler,
4851 : // doPointing, doPBCorr,
4852 : // tile_p, pbLimit_p, true);
4853 : // ((AWProjectFT *)ft_p)->setObservatoryLocation(mLocation_p);
4854 : // //
4855 : // // Explicit type casting of ft_p does not look good. It does not
4856 : // // pick up the setPAIncrement() method of PBWProjectFT without
4857 : // // this
4858 : // //
4859 : // Quantity paInc(paStep_p,"deg");
4860 : // // os << LogIO::NORMAL << "Setting PA increment to "
4861 : // // << paInc.getValue("deg") << " deg" << endl;
4862 : // ((AWProjectFT *)ft_p)->setPAIncrement(parAngleInc_p);
4863 :
4864 : // if (doPointing)
4865 : // {
4866 : // try
4867 : // {
4868 : // VisSet elVS(*rvi_p);
4869 : // epJ = new EPJones(elVS, *ms_p);
4870 : // RecordDesc applyRecDesc;
4871 : // applyRecDesc.addField("table", TpString);
4872 : // applyRecDesc.addField("interp",TpString);
4873 : // Record applyRec(applyRecDesc);
4874 : // applyRec.define("table",epJTableName_p);
4875 : // applyRec.define("interp", "nearest");
4876 : // epJ->setApply(applyRec);
4877 : // ((AWProjectFT *)ft_p)->setEPJones(epJ);
4878 : // }
4879 : // catch(AipsError& x)
4880 : // {
4881 : // //
4882 : // // Add some more useful info. to the message and translate
4883 : // // the generic AipsError exception object to a more specific
4884 : // // SynthesisError object.
4885 : // //
4886 : // String mesg = x.getMesg();
4887 : // mesg += ". Error in loading pointing offset table.";
4888 : // SynthesisError err(mesg);
4889 : // throw(err);
4890 : // }
4891 : // }
4892 : // AlwaysAssert(ft_p, AipsError);
4893 : // cft_p = new SimpleComponentFTMachine();
4894 : // AlwaysAssert(cft_p, AipsError);
4895 : // }
|
