Line data Source code
1 : //# SynthesisUtilMethods.cc:
2 : //# Copyright (C) 2013-2018
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 :
48 : #include <casacore/images/Images/TempImage.h>
49 : #include <casacore/images/Images/SubImage.h>
50 : #include <casacore/images/Regions/ImageRegion.h>
51 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
52 : #include <casacore/measures/Measures/MeasTable.h>
53 : #include <casacore/measures/Measures/MRadialVelocity.h>
54 : #include <casacore/ms/MSSel/MSSelection.h>
55 : #include <casacore/ms/MeasurementSets/MSColumns.h>
56 : #include <casacore/ms/MeasurementSets/MSDopplerUtil.h>
57 : #include <casacore/tables/Tables/Table.h>
58 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
59 : #include <synthesis/TransformMachines2/Utils.h>
60 :
61 : #include <mstransform/MSTransform/MSTransformRegridder.h>
62 : #include <msvis/MSVis/MSUtil.h>
63 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
64 : #include <msvis/MSVis/VisBufferUtil.h>
65 : #include <sys/types.h>
66 : #include <unistd.h>
67 : #include <limits>
68 : #include <tuple>
69 : #include <sys/time.h>
70 : #include<sys/resource.h>
71 :
72 : #include <synthesis/ImagerObjects/SIImageStore.h>
73 : #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
74 :
75 : using namespace std;
76 :
77 : using namespace casacore;
78 : namespace casa { //# NAMESPACE CASA - BEGIN
79 :
80 : casacore::String SynthesisUtilMethods::g_hostname;
81 : casacore::String SynthesisUtilMethods::g_startTimestamp;
82 : const casacore::String SynthesisUtilMethods::g_enableOptMemProfile =
83 : "synthesis.imager.memprofile.enable";
84 :
85 93 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 93 : }
89 :
90 93 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 93 : }
93 :
94 0 : Int SynthesisUtilMethods::validate(const VisBuffer& vb)
95 : {
96 0 : Int N=vb.nRow(),M=-1;
97 0 : for(Int i=0;i<N;i++)
98 : {
99 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
100 0 : {M++;break;}
101 : }
102 0 : return M;
103 : }
104 :
105 8382 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 8382 : Int N=vb.nRows(),M=-1;
108 26790 : for(Int i=0;i<N;i++)
109 : {
110 26790 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 8382 : {M++;break;}
112 : }
113 8382 : return M;
114 : }
115 : // Get the next largest even composite of 2,3,5.
116 : // This is to ensure a 'good' image size for FFTW.
117 : // Translated from gcwrap/scripts/cleanhelper.py : getOptimumSize
118 146 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 146 : Int n=npix;
121 :
122 146 : if( n%2 !=0 ){ n+= 1; }
123 :
124 146 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 1338 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 1192 : if( fac[k]>5 )
129 : {
130 0 : val = fac[k];
131 0 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 0 : fac[k] = val;
133 : }
134 : }
135 146 : newlarge=product(fac);
136 146 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 0 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 146 : return newlarge;
141 146 : }
142 :
143 : // Return the prime factors of the given number
144 146 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 146 : Vector<uInt> factors;
147 :
148 146 : Int lastresult = n;
149 146 : Int sqlast = int(sqrt(n))+1;
150 :
151 146 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 146 : Int c=2;
153 : while(1)
154 : {
155 1338 : if( lastresult == 1 || c > sqlast ) { break; }
156 1192 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 1192 : if( c>sqlast){ c=lastresult; break; }
160 1192 : if( lastresult % c == 0 ) { break; }
161 0 : c += 1;
162 : }
163 1192 : factors.resize( factors.nelements()+1, true );
164 1192 : factors[ factors.nelements()-1 ] = c;
165 1192 : lastresult /= c;
166 : }
167 146 : if( factors.nelements()==0 ) { factors.resize(1);factors[0]=n; }
168 :
169 : //if( douniq ) { factors = unique(factors); }
170 :
171 : /*
172 : /// The Sort::unique isn't working as called below. Results appear to be the
173 : /// same as with the cleanhelper python code, so leaving as is for not. CAS-7889
174 : if( douniq )
175 : {
176 : cout << "Test unique fn on : " << factors << endl;
177 : Sort srt;
178 : Vector<uInt> unvec=factors; uInt nrec;
179 : srt.unique(unvec,nrec);
180 : cout << " Unique : " << unvec << " nr : " << nrec << endl;
181 : }
182 : */
183 :
184 146 : return factors;
185 0 : }
186 :
187 :
188 0 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 0 : if (psfcutoff >=1.0 || psfcutoff<=0.0)
193 : {
194 0 : os << "psfcutoff must be >0 and <1" << LogIO::WARN;
195 0 : return false;
196 : }
197 :
198 0 : std::shared_ptr<SIImageStore> imstore;
199 0 : if( nterms>1 )
200 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 0 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 0 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 0 : imstore->printBeamSet();
210 :
211 0 : imstore->releaseLocks();
212 :
213 0 : return true;
214 0 : }
215 :
216 :
217 :
218 : // Evaluate a polynomial from Taylor coefficients and expand to a Cube.
219 : // Since this funcion is primarily for use with mtmfs_via_cube, this step applies only to
220 : // the multiterm.model.ttx images and the cube.model image cube.
221 0 : Bool SynthesisUtilMethods::taylorCoeffsToCube(const String& cubename,const String& mtname, const Int nterms, const String& reffreq)
222 : {
223 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "taylorCoeffsToCube"));
224 :
225 : // Set up imstores
226 0 : CountedPtr<SIImageStore> cube_imstore;
227 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
228 :
229 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
230 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
231 :
232 : // Check that .model exists.
233 : try{
234 0 : cube_imstore->model();
235 0 : for(Int i=0;i<nterms;i++)
236 0 : mt_imstore->model(i);
237 : }
238 0 : catch(AipsError &x)
239 : {
240 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\nModel images must exist on disk." ));
241 0 : }
242 :
243 : // Get/check shapes
244 0 : IPosition cube_shp( cube_imstore->model()->shape() );
245 0 : IPosition mt_shp( mt_imstore->model(0)->shape() );
246 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
247 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
248 : }
249 :
250 : // Read reference frequency
251 0 : Quantity reffreq_qa;
252 0 : Quantity::read( reffreq_qa, reffreq );
253 0 : Double refval = reffreq_qa.getValue("Hz");
254 :
255 : //cout << "ref freq : " << refval << endl;
256 :
257 : //Get the frequency list for the cube
258 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
259 0 : Vector<Double> freqlist( cube_shp[3] );
260 :
261 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
262 : {
263 0 : if( csys.type(i) == Coordinate::SPECTRAL )
264 : {
265 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
266 :
267 0 : for(Int ch=0;ch<cube_shp[3];ch++)
268 : {
269 : Double freq;
270 0 : Bool ret = speccoord.toWorld( freq, ch );
271 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
272 0 : freqlist[ch] = freq;
273 :
274 : //cout << "freq " << ch << " is " << freq << endl;
275 : }
276 :
277 0 : }
278 : }
279 :
280 :
281 : // Reset the Cube values to zero.
282 : //cube_imstore->model()->set(0.0);
283 :
284 : //For each pol, do the Taylor-to-Cube calculation.
285 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
286 : {
287 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(nterms);
288 0 : for(Int i=0;i<nterms;i++)
289 : {
290 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
291 0 : 0, cube_shp[3],
292 0 : pol, cube_shp[2],
293 0 : *mt_imstore->model(i) );
294 : }
295 :
296 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
297 : {
298 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
299 0 : chan, cube_shp[3],
300 0 : pol, cube_shp[2],
301 0 : *cube_imstore->model() );
302 :
303 0 : Double wt = (freqlist[chan] - refval) / refval;
304 :
305 0 : cube_subim->set(0.0);
306 0 : for(Int tt=0;tt<nterms;tt++)
307 : {
308 0 : Double fac = pow(wt,tt);
309 0 : LatticeExpr<Float> oneterm = LatticeExpr<Float>( *cube_subim + (fac) * (*mt_subims[tt])) ;
310 0 : cube_subim->copyData(oneterm);
311 0 : }
312 :
313 :
314 0 : }//for chan
315 0 : }// for pol
316 :
317 0 : return True;
318 :
319 0 : }//end of func
320 :
321 :
322 : // Calculate the RHS of the Normal equations for a linear least squares fit of a Taylor polynomial (per pixel).
323 : // This function is primarily for use with mtmfs_via_cube, and may be used for the PSF (2nterms-1), the residual (nterms)
324 : // and the primary beam (nterms=1).
325 : // imtype=0 : PSF with 2nterms-1 terms
326 : // imtype=1 : residual with nterms terms
327 : // imtype=2 : pb with 1 term
328 : // imtype=3 : sumwt with 2nterms-1 terms
329 0 : Bool SynthesisUtilMethods::cubeToTaylorSum(const String& cubename,const String& mtname, const Int nterms, const String& reffreq, const Int imtype, const Float pblimit)
330 : {
331 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "cubeToTaylorSum"));
332 :
333 : //cout << "imtype : " << imtype << endl;
334 0 : if(imtype <0 || imtype >3)
335 : {
336 0 : throw( AipsError("cubeToTaylorSum currently only supports 'psf','residual','pb', 'sumwt' options"));
337 : }
338 :
339 : // Set up imstores
340 0 : CountedPtr<SIImageStore> cube_imstore;
341 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
342 :
343 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
344 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
345 : //For psf avg pb has to be done already
346 : // doing residual too for sensitivity this is independent of beam spectral index removal
347 0 : Float maxPB=1.0;
348 0 : if(imtype < 2){
349 0 : LatticeExprNode elnod( max( *(mt_imstore->pb(0)) ) );
350 0 : maxPB=elnod.getFloat();
351 0 : if(maxPB == 0.0){
352 0 : throw(AipsError("Programmers error: should do tt psf images after making average PB"));
353 :
354 : }
355 :
356 0 : }
357 : // cerr << "imtype " << imtype << " MAX PB " << maxPB << endl;
358 : // If dopsf=True, calculate 2n-1 terms.
359 0 : Int out_nterms=nterms; // for residual
360 0 : if(imtype==0 || imtype==3){out_nterms=2 * nterms - 1;} // the psfs fill the upper triangle of the Hessian with 2 nterms-1 elements. Also sumwt.
361 0 : if(imtype==2){out_nterms=1;} // For the PB, for mtmfs_via_cube, we need only tt0. Later, if we need all terms to calculate PB alpha, then change this to nterms, and add the invHesian math (elsewhere) to later convert the RHS vector into the coefficients.
362 :
363 0 : CountedPtr <ImageInterface<Float> > use_cube, use_mt;
364 : // If dopsf=True, check that .psf cube and mt's exist. If dopsf=False, check residual images.
365 : try{
366 0 : switch(imtype)
367 : {
368 0 : case 0: use_cube=cube_imstore->psf();break;
369 0 : case 1: use_cube=cube_imstore->residual();break;
370 0 : case 2: use_cube=cube_imstore->pb();break;
371 0 : case 3: use_cube=cube_imstore->sumwt();break;
372 : }
373 0 : cube_imstore->sumwt();
374 0 : for(Int i=0;i<out_nterms;i++)
375 : {
376 0 : switch(imtype)
377 : {
378 0 : case 0:mt_imstore->psf(i);break;
379 0 : case 1:mt_imstore->residual(i);break;
380 0 : case 2:mt_imstore->pb(i);break;
381 0 : case 3:mt_imstore->sumwt(i);break;
382 : }
383 : }
384 : }
385 0 : catch(AipsError &x)
386 : {
387 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n " + imtype + " images must exist on disk." ));
388 0 : }
389 :
390 : // Get/check shapes ( Assume that the PSF always exists in the imstore... A valid assumption in the context of mtmfs_via_cube )
391 0 : IPosition cube_shp( cube_imstore->psf()->shape() );
392 0 : IPosition mt_shp( mt_imstore->psf(0)->shape() );
393 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
394 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
395 : }
396 :
397 : // Read reference frequency
398 0 : Quantity reffreq_qa;
399 0 : Quantity::read( reffreq_qa, reffreq );
400 0 : Double refval = reffreq_qa.getValue("Hz");
401 :
402 : //cout << "ref freq : " << refval << endl;
403 :
404 : //Get the frequency list for the cube
405 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
406 0 : Vector<Double> freqlist( cube_shp[3] );
407 :
408 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
409 : {
410 0 : if( csys.type(i) == Coordinate::SPECTRAL )
411 : {
412 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
413 :
414 0 : for(Int ch=0;ch<cube_shp[3];ch++)
415 : {
416 : Double freq;
417 0 : Bool ret = speccoord.toWorld( freq, ch );
418 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
419 0 : freqlist[ch] = freq;
420 : // cout << "freq " << ch << " is " << freq << endl;
421 : }
422 :
423 0 : }
424 : }
425 :
426 :
427 : // Reset the Taylor Sum values to zero.
428 0 : for(Int i=0;i<out_nterms;i++)
429 : {
430 0 : switch(imtype)
431 : {
432 0 : case 0:mt_imstore->psf(i)->set(0.0);break;
433 0 : case 1:mt_imstore->residual(i)->set(0.0);break;
434 0 : case 2:mt_imstore->pb(i)->set(0.0);break;
435 0 : case 3:mt_imstore->sumwt(i)->set(0.0);break;
436 : }
437 : }
438 :
439 : // Get the sumwt spectrum.
440 0 : Array<Float> lsumwt;
441 0 : cube_imstore->sumwt()->get(lsumwt, False);
442 :
443 : // Sum the weights ( or just use accumulate...)
444 0 : LatticeExprNode msum( sum( *cube_imstore->sumwt() ) );
445 0 : Float wtsum = msum.getFloat();
446 :
447 : //cerr << "perchansumwt : shape "<< lsumwt.shape() << " " << lsumwt << " sumwt "<< wtsum << endl;
448 :
449 : //Float wtsum = cube_shp[3]; // This is sum of weights, if all weights are 1.0
450 :
451 : //For each pol, do the Cube-To-Taylor calculation.
452 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
453 : {
454 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(out_nterms);
455 0 : for(Int i=0;i<out_nterms;i++)
456 : {
457 0 : switch(imtype)
458 : {
459 0 : case 0:use_mt=mt_imstore->psf(i);break;
460 0 : case 1:use_mt=mt_imstore->residual(i);break;
461 0 : case 2:use_mt=mt_imstore->pb(i);break;
462 0 : case 3:use_mt=mt_imstore->sumwt(i);break;
463 : }
464 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
465 0 : 0, cube_shp[3],
466 0 : pol, cube_shp[2],
467 0 : *use_mt );
468 : }
469 :
470 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
471 : {
472 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
473 0 : chan, cube_shp[3],
474 0 : pol, cube_shp[2],
475 0 : *use_cube);
476 0 : if(imtype < 2){
477 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
478 0 : chan, cube_shp[3],
479 0 : pol, cube_shp[2],
480 0 : *(cube_imstore->pb()));
481 0 : CountedPtr<ImageInterface<Float> > tmplat = new TempImage<Float>(cube_subim->shape(), cube_subim->coordinates());
482 0 : tmplat->copyData(LatticeExpr<Float>((*pb_subim) *(*cube_subim)));
483 0 : cube_subim = tmplat;
484 :
485 0 : }
486 :
487 0 : IPosition pos(4,0,0,pol,chan);
488 :
489 0 : Double wt = (freqlist[chan] - refval) / refval;
490 :
491 0 : for(Int tt=0;tt<out_nterms;tt++)
492 : {
493 0 : Double fac = pow(wt,tt);
494 : //cerr << "BEF accum " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
495 0 : LatticeExpr<Float> eachterm = LatticeExpr<Float>( (*mt_subims[tt]) + ((fac) * (*cube_subim) * lsumwt(pos))) ;
496 0 : mt_subims[tt]->copyData(eachterm);
497 : //cerr <<" AFT accum : chan " << chan << " tt " << tt << " fac " << fac << " lsumwt " << lsumwt(pos) << " pos " << pos << " max " << max(mt_subims[tt]->get()) << endl;
498 0 : }
499 :
500 :
501 0 : }//for chan
502 :
503 :
504 : // Divide by sum of weights.
505 0 : for(Int tt=0;tt<out_nterms;tt++)
506 : {
507 : //cerr << "bef div : tt " << tt << " : " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
508 :
509 0 : LatticeExpr<Float> eachterm;
510 0 : if (imtype < 2) {
511 0 : eachterm = LatticeExpr<Float>( iif( (*(mt_imstore->pb(0))) > pblimit , (*mt_subims[tt]) / wtsum/(*(mt_imstore->pb(0))), 0.0));
512 : }
513 : else{
514 0 : eachterm = LatticeExpr<Float>( (*mt_subims[tt]) / wtsum ) ;
515 : }
516 0 : mt_subims[tt]->copyData(eachterm);
517 :
518 0 : mt_subims[tt]->flush();
519 : // cerr << "aft div : " << max(mt_subims[tt]->get()) << endl;
520 0 : }
521 :
522 0 : }// for pol
523 :
524 :
525 : // Set the T/F mask, for PB images. Without this, the PB is fully masked, for aproj /mosaic gridders.
526 0 : if( imtype==2 )
527 : {
528 0 : mt_imstore->removeMask( mt_imstore->pb(0) );
529 : {
530 : //MSK//
531 0 : LatticeExpr<Bool> pbmask( iif( *mt_imstore->pb(0) > fabs(pblimit) , True , False ) );
532 : //MSK//
533 0 : mt_imstore->createMask( pbmask, mt_imstore->pb(0) );
534 0 : mt_imstore->pb(0)->pixelMask().unlock();
535 0 : }
536 :
537 : }
538 :
539 0 : return True;
540 :
541 0 : }//end of func
542 :
543 :
544 0 : Bool SynthesisUtilMethods::removeFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
545 : {
546 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "removeFreqDepPB"));
547 :
548 : // Set up imstores
549 0 : CountedPtr<SIImageStore> cube_imstore;
550 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
551 :
552 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
553 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
554 :
555 : try{
556 0 : cube_imstore->residual(); // Residual Cube
557 0 : cube_imstore->pb(); // PB cube
558 0 : mt_imstore->pb(0); // avgPB in the tt0 pb.
559 : }
560 0 : catch(AipsError &x)
561 : {
562 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n Residual cube, PB cube, and multiterm PB.tt0 must exist on disk." ));
563 0 : }
564 :
565 : // Get/check shapes
566 0 : IPosition cube_shp( cube_imstore->residual()->shape() );
567 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
568 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
569 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
570 : }
571 :
572 : //For each pol, do the freq-dep PB math.//////////////////////////
573 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
574 : {
575 :
576 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
577 0 : 0, cube_shp[3],
578 0 : pol, cube_shp[2],
579 0 : (*mt_imstore->pb(0)) );
580 :
581 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
582 0 : Float mtpbmaxval = mtpbmax.getFloat();
583 0 : if(mtpbmaxval <=0.0){os << LogIO::WARN << "pb.tt0 max is < or = zero" << LogIO::POST;}
584 :
585 :
586 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
587 : {
588 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
589 0 : chan, cube_shp[3],
590 0 : pol, cube_shp[2],
591 0 : *cube_imstore->residual() );
592 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
593 0 : chan, cube_shp[3],
594 0 : pol, cube_shp[2],
595 0 : *cube_imstore->pb() );
596 :
597 0 : LatticeExprNode pbmax( max( *pb_subim ) );
598 0 : Float pbmaxval = pbmax.getFloat();
599 0 : if( pbmaxval<=0.0 )
600 : {
601 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
602 : }
603 : else
604 : {
605 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*mt_subim)*(*(cube_subim))/(*(pb_subim)) , 0.0 ) );
606 0 : cube_subim->copyData( thepbcor );
607 0 : }// if not zero
608 :
609 0 : }//for chan
610 0 : }// for pol
611 :
612 0 : return True;
613 :
614 0 : }//end of func
615 :
616 :
617 :
618 0 : Bool SynthesisUtilMethods::applyFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
619 : {
620 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "applyFreqDepPB"));
621 :
622 : // Set up imstores
623 0 : CountedPtr<SIImageStore> cube_imstore;
624 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
625 :
626 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
627 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
628 :
629 : try{
630 0 : cube_imstore->model(); // Model Cube
631 0 : cube_imstore->pb(); // PB cube
632 0 : mt_imstore->pb(0); // avgPB in the tt0 pb.
633 : }
634 0 : catch(AipsError &x)
635 : {
636 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n Model cube, PB cube, and multiterm PB.tt0 must exist on disk." ));
637 0 : }
638 :
639 : // Get/check shapes
640 0 : IPosition cube_shp( cube_imstore->model()->shape() );
641 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
642 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
643 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
644 : }
645 :
646 : //For each pol, do the freq-dep PB math.//////////////////////////
647 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
648 : {
649 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
650 0 : 0, cube_shp[3],
651 0 : pol, cube_shp[2],
652 0 : (*mt_imstore->pb(0)) );
653 :
654 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
655 0 : Float mtpbmaxval = mtpbmax.getFloat();
656 0 : if(mtpbmaxval <=0.0)
657 0 : {os << LogIO::SEVERE << "pb.tt0 max is < or = zero. Cannot divide model image ! ERROR" << LogIO::POST;}
658 : else
659 : {
660 :
661 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
662 : {
663 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
664 0 : chan, cube_shp[3],
665 0 : pol, cube_shp[2],
666 0 : *cube_imstore->model() );
667 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
668 0 : chan, cube_shp[3],
669 0 : pol, cube_shp[2],
670 0 : *cube_imstore->pb() );
671 :
672 :
673 0 : LatticeExprNode pbmax( max( *pb_subim ) );
674 0 : Float pbmaxval = pbmax.getFloat();
675 0 : if( pbmaxval<=0.0 )
676 : {
677 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
678 : }
679 : else
680 : {
681 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*(cube_subim)) *(*(pb_subim)) / (*mt_subim) , 0.0 ) );
682 0 : cube_subim->copyData( thepbcor );
683 0 : }// if not zero
684 :
685 0 : }//for chan
686 : }// if mtpb >0
687 0 : }// for pol
688 :
689 0 : return True;
690 :
691 0 : }//end of func
692 :
693 :
694 :
695 :
696 : /***make a record of synthesisimager::weight parameters***/
697 73 : Record SynthesisUtilMethods::fillWeightRecord(const String& type, const String& rmode,
698 : const Quantity& noise, const Double robust,
699 : const Quantity& fieldofview,
700 : const Int npixels, const Bool multiField, const Bool useCubeBriggs,
701 : const String& filtertype, const Quantity& filterbmaj,
702 : const Quantity& filterbmin, const Quantity& filterbpa, const Double& fracBW){
703 :
704 73 : Record outRec;
705 73 : outRec.define("type", type);
706 73 : outRec.define("rmode", rmode);
707 73 : Record quantRec;
708 73 : QuantumHolder(noise).toRecord(quantRec);
709 73 : outRec.defineRecord("noise", quantRec);
710 73 : outRec.define("robust", robust);
711 73 : QuantumHolder(fieldofview).toRecord(quantRec);
712 73 : outRec.defineRecord("fieldofview", quantRec);
713 73 : outRec.define("npixels", npixels);
714 73 : outRec.define("multifield", multiField);
715 73 : outRec.define("usecubebriggs", useCubeBriggs);
716 73 : outRec.define("filtertype", filtertype);
717 73 : QuantumHolder(filterbmaj).toRecord(quantRec);
718 73 : outRec.defineRecord("filterbmaj", quantRec);
719 73 : QuantumHolder(filterbmin).toRecord(quantRec);
720 73 : outRec.defineRecord("filterbmin", quantRec);
721 73 : QuantumHolder(filterbpa).toRecord(quantRec);
722 73 : outRec.defineRecord("filterbpa", quantRec);
723 73 : outRec.define("fracBW", fracBW);
724 :
725 146 : return outRec;
726 73 : }
727 0 : void SynthesisUtilMethods::getFromWeightRecord(String& type, String& rmode,
728 : Quantity& noise, Double& robust,
729 : Quantity& fieldofview,
730 : Int& npixels, Bool& multiField, Bool& useCubeBriggs,
731 : String& filtertype, Quantity& filterbmaj,
732 : Quantity& filterbmin, Quantity& filterbpa, Double& fracBW, const Record& inRec){
733 0 : QuantumHolder qh;
734 0 : String err;
735 0 : if(!inRec.isDefined("type"))
736 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
737 0 : inRec.get("type", type);
738 0 : inRec.get("rmode", rmode);
739 0 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
740 0 : throw(AipsError("Error in reading noise param"));
741 0 : noise=qh.asQuantity();
742 0 : inRec.get("robust", robust);
743 0 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
744 0 : throw(AipsError("Error in reading fieldofview param"));
745 0 : fieldofview=qh.asQuantity();
746 0 : inRec.get("npixels", npixels);
747 0 : inRec.get("multifield", multiField);
748 0 : inRec.get("usecubebriggs", useCubeBriggs);
749 0 : inRec.get("filtertype", filtertype);
750 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
751 0 : throw(AipsError("Error in reading filterbmaj param"));
752 0 : filterbmaj=qh.asQuantity();
753 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
754 0 : throw(AipsError("Error in reading filterbmin param"));
755 0 : filterbmin=qh.asQuantity();
756 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
757 0 : throw(AipsError("Error in reading filterbpa param"));
758 0 : filterbpa=qh.asQuantity();
759 0 : inRec.get("fracBW", fracBW);
760 :
761 :
762 :
763 0 : }
764 :
765 :
766 : /**
767 : * Get values from lines of a /proc/self/status file. For example:
768 : * 'VmRSS: 600 kB'
769 : * @param str line from status file
770 : * @return integer value (memory amount, etc.)
771 : */
772 0 : Int SynthesisUtilMethods::parseProcStatusLine(const std::string &str) {
773 0 : istringstream is(str);
774 0 : std::string token;
775 0 : is >> token;
776 0 : is >> token;
777 0 : Int value = stoi(token);
778 0 : return value;
779 0 : }
780 :
781 : /**
782 : * Produces a name for a 'memprofile' output file. For example:
783 : * casa.synthesis.imager.memprofile.23514.pc22555.hq.eso.org.20171209_120446.txt
784 : * (where 23514 is the PID passed as input parameter).
785 : *
786 : * @param pid PID of the process running the imager
787 : *
788 : * @return a longish 'memprofile' filename including PID, machine, timestamp, etc.
789 : **/
790 0 : String SynthesisUtilMethods::makeResourceFilename(int pid)
791 : {
792 0 : if (g_hostname.empty() or g_startTimestamp.empty()) {
793 : // TODO: not using HOST_NAME_MAX because of issues with __APPLE__
794 : // somehow tests tAWPFTM, tSynthesisImager, and tSynthesisImagerVi2 fail.
795 0 : const int strMax = 255;
796 : char hostname[strMax];
797 0 : gethostname(hostname, strMax);
798 0 : g_hostname = hostname;
799 :
800 0 : auto time = std::time(nullptr);
801 0 : auto gmt = std::gmtime(&time);
802 0 : const char* format = "%Y%m%d_%H%M%S";
803 : char timestr[strMax];
804 0 : std::strftime(timestr, strMax, format, gmt);
805 0 : g_startTimestamp = timestr;
806 : }
807 :
808 0 : return String("casa.synthesis.imager.memprofile." + String::toString(pid) +
809 0 : "." + g_hostname + "." + g_startTimestamp + ".txt");
810 : }
811 :
812 0 : Bool SynthesisUtilMethods::adviseChanSel(Double& freqStart, Double& freqEnd,
813 : const Double& freqStep, const MFrequency::Types& freqframe,
814 : Vector<Int>& spw, Vector<Int>& start,
815 : Vector<Int>& nchan, const String& ms, const String& ephemtab, const Int field_id, const Bool getFreqRange, const String spwselection){
816 :
817 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "adviseChanSel"));
818 0 : if(ms==String("")){
819 0 : throw(AipsError("Need a valid MS"));
820 : }
821 0 : spw.resize();
822 0 : start.resize();
823 0 : nchan.resize();
824 : try {
825 0 : if(!getFreqRange){
826 0 : Vector<Int> bnchan;
827 0 : Vector<Int> bstart;
828 0 : Vector<Int> bspw;
829 : Double fS, fE;
830 0 : fS=freqStart;
831 0 : fE=freqEnd;
832 0 : if(freqEnd < freqStart){
833 0 : fS=freqEnd;
834 0 : fE=freqStart;
835 : }
836 :
837 :
838 : {
839 :
840 0 : MeasurementSet elms(String(ms), TableLock(TableLock::AutoNoReadLocking), Table::Old);
841 0 : if(ephemtab != "" && freqframe == MFrequency::REST ){
842 0 : MSUtil::getSpwInSourceFreqRange(bspw, bstart, bnchan, elms, fS, fE, fabs(freqStep), ephemtab, field_id);
843 : }
844 : else
845 0 : MSUtil::getSpwInFreqRange(bspw, bstart, bnchan, elms, fS, fE, fabs(freqStep), freqframe, field_id);
846 0 : elms.relinquishAutoLocks(true);
847 :
848 0 : }
849 0 : spw=Vector<Int> (bspw);
850 0 : start=Vector<Int> (bstart);
851 0 : nchan=Vector<Int> (bnchan);
852 0 : }
853 : else{
854 :
855 : {
856 0 : MeasurementSet elms(ms, TableLock(TableLock::AutoNoReadLocking), Table::Old);
857 0 : MSSelection thisSelection;
858 0 : String spsel=spwselection;
859 0 : if(spsel=="")spsel="*";
860 0 : thisSelection.setSpwExpr(spsel);
861 0 : TableExprNode exprNode=thisSelection.toTableExprNode(&elms);
862 0 : Matrix<Int> chanlist=thisSelection.getChanList();
863 0 : if(chanlist.ncolumn() <3){
864 0 : freqStart=-1.0;
865 0 : freqEnd=-1.0;
866 0 : return false;
867 : }
868 0 : Vector<Int> elspw=chanlist.column(0);
869 0 : Vector<Int> elstart=chanlist.column(1);
870 0 : Vector<Int> elnchan=Vector<Int> (chanlist.column(2)-elstart)+1;
871 0 : if(ephemtab != "" ){
872 0 : const MSColumns mscol(ms);
873 0 : MEpoch ep=mscol.timeMeas()(0);
874 0 : Quantity sysvel;
875 0 : String ephemTable("");
876 0 : MDirection::Types mtype=MDirection::APP;
877 0 : MDirection mdir(mtype);
878 0 : if(Table::isReadable(ephemtab)){
879 0 : ephemTable=ephemtab;
880 : }
881 0 : else if(ephemtab=="TRACKFIELD"){
882 0 : ephemTable=(mscol.field()).ephemPath(field_id);
883 : }
884 0 : else if(MDirection::getType(mtype, ephemtab)){
885 0 : mdir=MDirection(mtype);
886 : }
887 :
888 0 : MSUtil::getFreqRangeAndRefFreqShift(freqStart, freqEnd, sysvel, ep, elspw, elstart, elnchan, elms, ephemTable , mdir, True);
889 :
890 0 : }
891 : else
892 0 : MSUtil::getFreqRangeInSpw(freqStart, freqEnd, elspw, elstart, elnchan, elms, freqframe, field_id);
893 0 : }
894 :
895 : }
896 :
897 :
898 :
899 :
900 0 : } catch (AipsError x) {
901 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
902 0 : << LogIO::POST;
903 0 : return false;
904 0 : }
905 0 : catch (...){
906 : os << LogIO::SEVERE << "Unknown exception handled"
907 0 : << LogIO::POST;
908 0 : return false;
909 :
910 0 : }
911 :
912 0 : return true;
913 :
914 0 : }
915 :
916 3152 : void SynthesisUtilMethods::getResource(String label, String fname)
917 : {
918 : // TODO: not tested on anything else than LINUX (MACOS for the future)
919 : #if !defined(AIPS_LINUX)
920 : return;
921 : #endif
922 :
923 3152 : Bool isOn = false;
924 3152 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
925 3152 : if (!isOn)
926 3152 : return;
927 :
928 : // TODO: reorganize, use struct or something to hold and pass info over. ifdef lnx
929 0 : LogIO casalog( LogOrigin("SynthesisUtilMethods", "getResource", WHERE) );
930 :
931 : // To hold memory stats, in MB
932 0 : int vmRSS = -1, vmHWM = -1, vmSize = -1, vmPeak = -1, vmSwap = -1;
933 0 : pid_t pid = -1;
934 0 : int fdSize = -1;
935 :
936 : // TODO: this won't probably work on anything but linux
937 0 : ifstream procFile("/proc/self/status");
938 0 : if (procFile.is_open()) {
939 0 : std::string line;
940 0 : while (not procFile.eof()) {
941 0 : getline(procFile, line);
942 0 : const std::string startVmRSS = "VmRSS:";
943 0 : const std::string startVmWHM = "VmHWM:";
944 0 : const std::string startVmSize = "VmSize:";
945 0 : const std::string startVmPeak = "VmPeak:";
946 0 : const std::string startVmSwap = "VmSwap:";
947 0 : const std::string startFDSize = "FDSize:";
948 0 : const double KB_TO_MB = 1024.0;
949 0 : if (startVmRSS == line.substr(0, startVmRSS.size())) {
950 0 : vmRSS = parseProcStatusLine(line.c_str()) / KB_TO_MB;
951 0 : } else if (startVmWHM == line.substr(0, startVmWHM.size())) {
952 0 : vmHWM = parseProcStatusLine(line.c_str()) / KB_TO_MB;
953 0 : } else if (startVmSize == line.substr(0, startVmSize.size())) {
954 0 : vmSize = parseProcStatusLine(line.c_str()) / KB_TO_MB;
955 0 : } else if (startVmPeak == line.substr(0, startVmPeak.size())) {
956 0 : vmPeak = parseProcStatusLine(line.c_str()) / KB_TO_MB;
957 0 : } else if (startVmSwap == line.substr(0, startVmSwap.size())) {
958 0 : vmSwap = parseProcStatusLine(line.c_str()) / KB_TO_MB;
959 0 : } else if (startFDSize == line.substr(0, startFDSize.size())) {
960 0 : fdSize = parseProcStatusLine(line.c_str());
961 : }
962 0 : }
963 0 : procFile.close();
964 0 : }
965 :
966 0 : pid = getpid();
967 :
968 : struct rusage usage;
969 : struct timeval now;
970 0 : getrusage(RUSAGE_SELF, &usage);
971 0 : now = usage.ru_utime;
972 :
973 : // TODO: check if this works as expected when /proc/self/status is not there
974 : // Not clear at all if VmHWM and .ru_maxrss measure the same thing
975 : // Some alternative is needed for the other fields as well: VmSize, VMHWM, FDSize.
976 0 : if (vmHWM < 0) {
977 0 : vmHWM = usage.ru_maxrss;
978 : }
979 :
980 0 : ostringstream oss;
981 0 : oss << "PID: " << pid ;
982 0 : oss << " MemRSS (VmRSS): " << vmRSS << " MB.";
983 0 : oss << " VmWHM: " << vmHWM << " MB.";
984 0 : oss << " VirtMem (VmSize): " << vmSize << " MB.";
985 0 : oss << " VmPeak: " << vmPeak << " MB.";
986 0 : oss << " VmSwap: " << vmSwap << " MB.";
987 0 : oss << " ProcTime: " << now.tv_sec << '.' << now.tv_usec;
988 0 : oss << " FDSize: " << fdSize;
989 0 : oss << " [" << label << "] ";
990 0 : casalog << oss.str() << LogIO::NORMAL3 << LogIO::POST;
991 :
992 : // Write this to a file too...
993 : try {
994 0 : if (fname.empty()) {
995 0 : fname = makeResourceFilename(pid);
996 : }
997 0 : ofstream ofile(fname, ios::app);
998 0 : if (ofile.is_open()) {
999 0 : if (0 == ofile.tellp()) {
1000 : casalog << g_enableOptMemProfile << " is enabled, initializing output file for "
1001 : "imager profiling information (memory and run time): " << fname <<
1002 0 : LogIO::NORMAL << LogIO::POST;
1003 0 : ostringstream header;
1004 : header << "# PID, MemRSS_(VmRSS)_MB, VmWHM_MB, VirtMem_(VmSize)_MB, VmPeak_MB, "
1005 0 : "VmSwap_MB, ProcTime_sec, FDSize, label_checkpoint";
1006 0 : ofile << header.str() << '\n';
1007 0 : }
1008 0 : ostringstream line;
1009 0 : line << pid << ',' << vmRSS << ',' << vmHWM << ',' << vmSize << ','
1010 0 : << vmPeak << ','<< vmSwap << ',' << now.tv_sec << '.' << now.tv_usec << ','
1011 0 : << fdSize << ',' << '[' << label << ']';
1012 0 : ofile << line.str() << '\n';
1013 0 : ofile.close();
1014 0 : }
1015 0 : } catch(std::runtime_error &exc) {
1016 : casalog << "Could not write imager memory+runtime information into output file: "
1017 0 : << fname << LogIO::WARN << LogIO::POST;
1018 0 : }
1019 0 : }
1020 :
1021 :
1022 : // Data partitioning rules for CONTINUUM imaging
1023 : //
1024 : // ALL members of the selection parameters in selpars are strings
1025 : // ONLY. This methods reads the selection parameters from selpars
1026 : // and returns a partitioned Record with npart data selection
1027 : // entries.
1028 : //
1029 : // The algorithm used to do the partitioning along the TIME axis is
1030 : // as follows:
1031 : //
1032 : // for each MS in selpars
1033 : // - get the selection parameters
1034 : // - generate a selected MS
1035 : // - get number of rows in the selected MS
1036 : // - divide the rows in nparts
1037 : // - for each part
1038 : // - get compute rowBeg and rowEnd
1039 : // - modify rowEnd such that rowEnd points to the end of
1040 : // full integration data. This is done as follows:
1041 : // tRef = TIME(rowEnd);
1042 : // reduce rowEnd till TIME(rowEnd) != tRef
1043 : // - Construct a T0~T1 string
1044 : // - Fill it in the timeSelPerPart[msID][PartNo] array
1045 : //
1046 0 : Record SynthesisUtilMethods::continuumDataPartition(Record &selpars, const Int npart)
1047 : {
1048 0 : LogIO os( LogOrigin("SynthesisUtilMethods","continuumDataPartition",WHERE) );
1049 :
1050 0 : Record onepart, allparts;
1051 0 : Vector<Vector<String> > timeSelPerPart;
1052 0 : timeSelPerPart.resize(selpars.nfields());
1053 :
1054 : // Duplicate the entire input record npart times, with a specific partition id.
1055 : // Modify each sub-record according to partition id.
1056 0 : for (uInt msID=0;msID<selpars.nfields();msID++)
1057 : {
1058 0 : Record thisMS= selpars.subRecord(RecordFieldId("ms"+String::toString(msID)));
1059 0 : String msName = thisMS.asString("msname");
1060 0 : timeSelPerPart[msID].resize(npart,true);
1061 : //
1062 : // Make a selected MS and extract the time-column information
1063 : //
1064 0 : MeasurementSet ms(msName,TableLock(TableLock::AutoNoReadLocking), Table::Old),
1065 0 : selectedMS(ms);
1066 0 : MSInterface msI(ms); MSSelection msSelObj;
1067 0 : msSelObj.reset(msI,MSSelection::PARSE_NOW,
1068 : thisMS.asString("timestr"),
1069 : thisMS.asString("antenna"),
1070 : thisMS.asString("field"),
1071 : thisMS.asString("spw"),
1072 : thisMS.asString("uvdist"),
1073 : thisMS.asString("taql"),
1074 : "",//thisMS.asString("poln"),
1075 : thisMS.asString("scan"),
1076 : "",//thisMS.asString("array"),
1077 : thisMS.asString("state"),
1078 : thisMS.asString("obs")//observation
1079 : );
1080 0 : msSelObj.getSelectedMS(selectedMS);
1081 :
1082 : //--------------------------------------------------------------------
1083 : // Use the selectedMS to generate time selection strings per part
1084 : //
1085 : // Double Tint;
1086 0 : MSMainColumns mainCols(selectedMS);
1087 0 : Vector<rownr_t> rowNumbers = selectedMS.rowNumbers();
1088 0 : Int nRows=selectedMS.nrow(),
1089 0 : dRows=nRows/npart;
1090 0 : Int rowBegID=0, rowEndID=nRows-1;
1091 0 : Int rowBeg=rowNumbers[rowBegID], rowEnd=rowNumbers[rowEndID];
1092 : //cerr << "NRows, dRows, npart = " << nRows << " " << dRows << " " << npart << " " << rowBeg << " " << rowEnd << endl;
1093 :
1094 0 : rowEndID = rowBegID + dRows;
1095 :
1096 :
1097 0 : MVTime mvInt=mainCols.intervalQuant()(0);
1098 : //Time intT(mvInt.getTime());
1099 : // Tint = intT.modifiedJulianDay();
1100 :
1101 0 : Int partNo=0;
1102 : // The +1 in rowBeg and rowEnd below is because it appears
1103 : // that TaQL by default counts from 1, not 0.
1104 0 : while(rowEndID < nRows)
1105 : {
1106 : // rowBeg=rowNumbers[rowBegID]; rowEnd = rowNumbers[rowEndID];
1107 0 : rowBeg=rowBegID+1; rowEnd = rowEndID+1;
1108 0 : stringstream taql;
1109 0 : taql << "ROWNUMBER() >= " << rowBeg << " && ROWNUMBER() <= " << rowEnd;
1110 0 : timeSelPerPart[msID][partNo] = taql.str();
1111 :
1112 0 : if (partNo == npart - 1) break;
1113 0 : partNo++;
1114 0 : rowBegID = rowEndID+1;
1115 0 : rowEndID = min(rowBegID + dRows, nRows-1);
1116 0 : if (rowEndID == nRows-1) break;
1117 0 : }
1118 :
1119 : //rowBeg=rowNumbers[rowBegID]; rowEnd = rowNumbers[nRows-1];
1120 0 : stringstream taql;
1121 0 : rowBeg=rowBegID+1; rowEnd = nRows-1+1;
1122 0 : taql << "ROWNUMBER() >= " << rowBeg << " && ROWNUMBER() <= " << rowEnd;
1123 0 : timeSelPerPart[msID][partNo] = taql.str();
1124 : os << endl << "Rows = " << rowBeg << " " << rowEnd << " "
1125 0 : << "[P][M]: " << msID << ":" << partNo << " " << timeSelPerPart[msID][partNo]
1126 0 : << LogIO::POST;
1127 0 : }
1128 : //
1129 : // The time selection strings for all parts of the current
1130 : // msID are in timeSelPerPart.
1131 : //--------------------------------------------------------------------
1132 : //
1133 : // Now reverse the order of parts and ME loops. Create a Record
1134 : // per part, get the MS for thisPart. Put the associated time
1135 : // selection string in it. Add the thisMS to thisPart Record.
1136 : // Finally add thisPart Record to the allparts Record.
1137 : //
1138 0 : for(Int iPart=0; iPart<npart; iPart++)
1139 : {
1140 0 : Record thisPart;
1141 0 : thisPart.assign(selpars);
1142 0 : for (uInt msID=0; msID<selpars.nfields(); msID++)
1143 : {
1144 0 : Record thisMS = thisPart.subRecord(RecordFieldId("ms"+String::toString(msID)));
1145 :
1146 0 : thisMS.define("taql",timeSelPerPart[msID][iPart]);
1147 0 : thisPart.defineRecord(RecordFieldId("ms"+String::toString(msID)) , thisMS);
1148 0 : }
1149 0 : allparts.defineRecord(RecordFieldId(String::toString(iPart)), thisPart);
1150 0 : }
1151 : // cerr << allparts << endl;
1152 0 : return allparts;
1153 :
1154 : // for( Int part=0; part < npart; part++)
1155 : // {
1156 :
1157 : // onepart.assign(selpars);
1158 :
1159 :
1160 : // //-------------------------------------------------
1161 : // // WARNING : This is special-case code for two specific datasets
1162 : // for ( uInt msid=0; msid<selpars.nfields(); msid++)
1163 : // {
1164 : // Record onems = onepart.subRecord( RecordFieldId("ms"+String::toString(msid)) );
1165 : // String msname = onems.asString("msname");
1166 : // String spw = onems.asString("spw");
1167 : // if(msname.matches("DataTest/twopoints_twochan.ms"))
1168 : // {
1169 : // onems.define("spw", spw+":"+String::toString(part));
1170 : // }
1171 : // if(msname.matches("DataTest/point_twospws.ms"))
1172 : // {
1173 : // onems.define("spw", spw+":"+ (((Bool)part)?("10~19"):("0~9")) );
1174 : // }
1175 : // if(msname.matches("DataTest/reg_mawproject.ms"))
1176 : // {
1177 : // onems.define("scan", (((Bool)part)?("1~17"):("18~35")) );
1178 : // }
1179 : // onepart.defineRecord( RecordFieldId("ms"+String::toString(msid)) , onems );
1180 : // }// end ms loop
1181 : // //-------------------------------------------------
1182 :
1183 : // allparts.defineRecord( RecordFieldId(String::toString(part)), onepart );
1184 :
1185 : // }// end partition loop
1186 :
1187 : // return allparts;
1188 0 : }
1189 :
1190 :
1191 : // Data partitioning rules for CUBE imaging
1192 0 : Record SynthesisUtilMethods::cubeDataPartition(const Record &selpars, const Int npart,
1193 : const Double freqBeg, const Double freqEnd, const MFrequency::Types eltype)
1194 : {
1195 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataPartition",WHERE) );
1196 : // Temporary special-case code. Please replace with actual rules.
1197 0 : Vector<Double> fstart(npart);
1198 0 : Vector<Double> fend(npart);
1199 0 : Double step=(freqEnd-freqBeg)/Double(npart);
1200 0 : fstart(0)=freqBeg;
1201 0 : fend(0)=freqBeg+step;
1202 0 : for (Int k=1; k < npart; ++k){
1203 0 : fstart(k)=fstart(k-1)+step;
1204 0 : fend(k)=fend(k-1)+step;
1205 : }
1206 0 : return cubeDataPartition( selpars, fstart, fend, eltype );
1207 :
1208 0 : }
1209 :
1210 :
1211 0 : Record SynthesisUtilMethods::cubeDataImagePartition(const Record & selpars, const CoordinateSystem&
1212 : incsys, const Int npart, const Int nchannel,
1213 : Vector<CoordinateSystem>& outCsys,
1214 : Vector<Int>& outnChan){
1215 :
1216 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataImagePartition",WHERE) );
1217 0 : outnChan.resize(npart);
1218 0 : outCsys.resize(npart);
1219 0 : Int nomnchan=nchannel/npart;
1220 0 : outnChan.set(nomnchan);
1221 0 : nomnchan=nchannel%npart;
1222 0 : for (Int k=0; k < nomnchan; ++k)
1223 0 : outnChan[k]+=1;
1224 0 : Vector<Int> shp(0);
1225 : //shp(0)=20; shp(1)=20; shp(2)=1; shp(3)=outnChan[0];
1226 0 : Vector<Float> shift(4, 0.0);
1227 0 : Vector<Float> fac(4, 1.0);
1228 0 : Vector<Double> freqEnd(npart);
1229 0 : Vector<Double> freqStart(npart);
1230 0 : Float chanshift=0.0;
1231 0 : for (Int k =0; k <npart; ++k){
1232 0 : shift(3)=chanshift;
1233 : //cerr << k << " shift " << shift << endl;
1234 0 : outCsys[k]=incsys.subImage(shift, fac, shp);
1235 0 : freqStart[k]=SpectralImageUtil::worldFreq(outCsys[k], 0.0);
1236 0 : freqEnd[k]=SpectralImageUtil::worldFreq(outCsys[k], Double(outnChan[k]-1));
1237 0 : if(freqStart[k] > freqEnd[k]){
1238 0 : Double tmp=freqEnd[k];
1239 0 : freqEnd[k]=freqStart[k];
1240 0 : freqStart[k]=tmp;
1241 : }
1242 0 : chanshift+=Float(outnChan[k]);
1243 : }
1244 0 : MFrequency::Types eltype=incsys.spectralCoordinate(incsys.findCoordinate(Coordinate::SPECTRAL)).frequencySystem(true);
1245 :
1246 : //os << "freqStart="<<freqStart<<" freqend="<<freqEnd<< "eltype="<<eltype<<LogIO::POST;
1247 0 : Record rec=cubeDataPartition(selpars, freqStart, freqEnd, eltype);
1248 0 : for (Int k=0; k < npart ; ++k){
1249 0 : outCsys[k].save(rec.asrwRecord(String::toString(k)), "coordsys");
1250 0 : rec.asrwRecord(String::toString(k)).define("nchan", outnChan[k]);
1251 : }
1252 0 : return rec;
1253 0 : }
1254 :
1255 0 : Record SynthesisUtilMethods::cubeDataPartition(const Record &selpars, const Vector<Double>& freqBeg, const Vector<Double>&freqEnd, const MFrequency::Types eltype){
1256 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataPartition",WHERE) );
1257 0 : Record retRec;
1258 0 : Int npart=freqBeg.shape()(0);
1259 0 : for (Int k=0; k < npart; ++k){
1260 0 : Int nms=selpars.nfields();
1261 0 : Record partRec;
1262 0 : for(Int j=0; j < nms; ++j){
1263 0 : if(selpars.isDefined(String("ms"+String::toString(j)))){
1264 0 : Record msRec=selpars.asRecord(String("ms"+String::toString(j)));
1265 0 : if(!msRec.isDefined("msname"))
1266 0 : throw(AipsError("No msname key in ms record"));
1267 0 : String msname=msRec.asString("msname");
1268 0 : String userspw=msRec.isDefined("spw")? msRec.asString("spw") : "*";
1269 0 : String userfield=msRec.isDefined("field") ? msRec.asString("field") : "*";
1270 0 : String userstate=msRec.isDefined("state") ? msRec.asString("state") : "*";
1271 :
1272 0 : MeasurementSet elms(msname);
1273 0 : Record laSelection=elms.msseltoindex(userspw, userfield);
1274 0 : if (userfield=="") userfield="*";
1275 0 : MSSelection mssel;
1276 0 : mssel.setSpwExpr(userspw);
1277 0 : mssel.setFieldExpr(userfield);
1278 0 : mssel.setStateExpr(userstate);
1279 0 : TableExprNode exprNode = mssel.toTableExprNode(&elms);
1280 0 : Matrix<Int> spwsel=mssel.getChanList();
1281 0 : Vector<Int> fieldsel=mssel.getFieldList();
1282 : // case for scan intent specified
1283 0 : if (userstate!="*") {
1284 0 : MeasurementSet elselms((elms)(exprNode), &elms);
1285 0 : MSColumns tmpmsc(elselms);
1286 0 : Vector<Int> fldidv=tmpmsc.fieldId().getColumn();
1287 0 : if (fldidv.nelements()==0)
1288 0 : throw(AipsError("No field ids were selected, please check input parameters"));
1289 0 : std::set<Int> ufldids(fldidv.begin(),fldidv.end());
1290 0 : std::vector<Int> tmpv(ufldids.begin(), ufldids.end());
1291 0 : fieldsel.resize(tmpv.size());
1292 0 : uInt count=0;
1293 0 : for (std::vector<int>::const_iterator it=tmpv.begin();it != tmpv.end(); it++)
1294 : {
1295 0 : fieldsel(count) = *it;
1296 0 : count++;
1297 : }
1298 0 : }
1299 : //Matrix<Int> spwsel=laSelection.asArrayInt("channel");
1300 : //Vector<Int> fieldsel=laSelection.asArrayInt("field");
1301 0 : Vector<Int> freqSpw;
1302 0 : Vector<Int> freqStart;
1303 0 : Vector<Int> freqNchan;
1304 0 : String newspw;
1305 : try{
1306 0 : MSUtil::getSpwInFreqRange(freqSpw, freqStart, freqNchan, elms, freqBeg(k), freqEnd(k),0.0, eltype, fieldsel[0]);
1307 0 : newspw=mergeSpwSel(freqSpw, freqStart, freqNchan, spwsel);
1308 : //cerr << "try " << freqSpw << " " << freqStart << " " << freqNchan << endl;
1309 : }
1310 0 : catch(...){
1311 : //cerr << "In catch " << endl;
1312 0 : newspw="";
1313 0 : }
1314 : //String newspw=mergeSpwSel(freqSpw, freqStart, freqNchan, spwsel);
1315 0 : if(newspw=="") newspw="-1";
1316 0 : msRec.define("spw", newspw);
1317 0 : partRec.defineRecord(String("ms"+String::toString(j)),msRec);
1318 0 : }
1319 :
1320 : }
1321 0 : retRec.defineRecord(String::toString(k), partRec);
1322 0 : }
1323 :
1324 :
1325 :
1326 :
1327 0 : return retRec;
1328 0 : }
1329 :
1330 :
1331 0 : String SynthesisUtilMethods::mergeSpwSel(const Vector<Int>& fspw, const Vector<Int>& fstart, const Vector<Int>& fnchan, const Matrix<Int>& spwsel)
1332 : {
1333 0 : String retval="";
1334 : Int cstart, cend;
1335 0 : for(Int k=0; k < fspw.shape()(0); ++k){
1336 0 : cstart=fstart(k);
1337 0 : cend=fstart(k)+fnchan(k)-1;
1338 0 : for (Int j=0; j < spwsel.shape()(0); ++j){
1339 : //need to put this here as multiple rows can have the same spw
1340 0 : cstart=fstart(k);
1341 0 : cend=fstart(k)+fnchan(k)-1;
1342 0 : if(spwsel(j,0)==fspw[k]){
1343 0 : if(cstart < spwsel(j,1)) cstart=spwsel(j,1);
1344 0 : if(cend > spwsel(j,2)) cend= spwsel(j,2);
1345 0 : if(!retval.empty()) retval=retval+(",");
1346 0 : retval=retval+String::toString(fspw[k])+":"+String::toString(cstart)+"~"+String::toString(cend);
1347 : }
1348 : }
1349 : }
1350 :
1351 :
1352 :
1353 0 : return retval;
1354 0 : }
1355 :
1356 : // Image cube partitioning rules for CUBE imaging
1357 0 : Record SynthesisUtilMethods::cubeImagePartition(Record &impars, Int npart)
1358 : {
1359 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeImagePartition",WHERE) );
1360 :
1361 0 : Record onepart, allparts;
1362 :
1363 : // Duplicate the entire input record npart times, with a specific partition id.
1364 : // Modify each sub-record according to partition id.
1365 0 : for( Int part=0; part < npart; part++)
1366 : {
1367 :
1368 0 : onepart.assign(impars);
1369 :
1370 : //-------------------------------------------------
1371 : // WARNING : This is special-case code for two specific datasets
1372 0 : for ( uInt imfld=0; imfld<impars.nfields(); imfld++)
1373 : {
1374 0 : Record onefld = onepart.subRecord( RecordFieldId(String::toString(imfld)) );
1375 0 : Int nchan = onefld.asInt("nchan");
1376 : //String freqstart = onems.asString("freqstart");
1377 :
1378 0 : onefld.define("nchan", nchan/npart);
1379 0 : onefld.define("freqstart", (((Bool)part)?("1.2GHz"):("1.0GHz")) );
1380 :
1381 0 : String imname = onefld.asString("imagename");
1382 0 : onefld.define("imagename", imname+".n"+String::toString(part));
1383 :
1384 0 : onepart.defineRecord( RecordFieldId( String::toString(imfld) ), onefld );
1385 0 : }// end ms loop
1386 : //-------------------------------------------------
1387 :
1388 0 : allparts.defineRecord( RecordFieldId(String::toString(part)), onepart );
1389 :
1390 : }// end partition loop
1391 :
1392 0 : return allparts;
1393 :
1394 :
1395 0 : }
1396 :
1397 0 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
1398 : {
1399 0 : MVAngle mvRa=direction.getAngle().getValue()(0);
1400 0 : MVAngle mvDec=direction.getAngle().getValue()(1);
1401 0 : ostringstream oos;
1402 0 : oos << " ";
1403 0 : Int widthRA=20;
1404 0 : Int widthDec=20;
1405 0 : oos.setf(ios::left, ios::adjustfield);
1406 0 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
1407 0 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
1408 0 : oos << " "
1409 0 : << MDirection::showType(direction.getRefPtr()->getType());
1410 0 : return String(oos);
1411 0 : }
1412 :
1413 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1414 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1415 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
1416 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1417 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1418 :
1419 : // Read String from Record
1420 4745 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1421 : {
1422 4745 : if( rec.isDefined( id ) )
1423 : {
1424 4380 : String inval("");
1425 4380 : if( rec.dataType( id )==TpString )
1426 4380 : { rec.get( id , inval ); // Read into temp string
1427 : // val = inval;
1428 : // return String("");
1429 : // Set value only if it is not a null string. Otherwise, leave it unchanged as it will
1430 : // retain the default value that was set before this function was called.
1431 4380 : if(inval.length()>0){val=inval;}
1432 4380 : return String("");
1433 : }
1434 0 : else { return String(id + " must be a string\n"); }
1435 4380 : }
1436 365 : else { return String("");}
1437 : }
1438 :
1439 : // Read Integer from Record
1440 1170 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1441 : {
1442 1170 : if( rec.isDefined( id ) )
1443 : {
1444 1168 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1445 0 : else { return String(id + " must be an integer\n"); }
1446 : }
1447 2 : else { return String(""); }
1448 : }
1449 :
1450 : // Read Float from Record
1451 2117 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1452 : {
1453 2117 : if( rec.isDefined( id ) )
1454 : {
1455 1971 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1456 1971 : { rec.get( id , val ); return String(""); }
1457 0 : else { return String(id + " must be a float\n"); }
1458 : }
1459 146 : else { return String(""); }
1460 : }
1461 :
1462 : // Read Bool from Record
1463 2263 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1464 : {
1465 2263 : if( rec.isDefined( id ) )
1466 : {
1467 1752 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1468 0 : else { return String(id + " must be a bool\n"); }
1469 : }
1470 511 : else{ return String(""); }
1471 : }
1472 :
1473 : // Read Vector<Int> from Record
1474 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
1475 : {
1476 0 : if( rec.isDefined( id ) )
1477 : {
1478 0 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1479 0 : { rec.get( id , val ); return String(""); }
1480 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1481 : {
1482 0 : Vector<Bool> vec; rec.get( id, vec );
1483 0 : if( vec.nelements()==0 ){val.resize(0); return String("");}
1484 0 : else{ return String(id + " must be a vector of strings.\n"); }
1485 0 : }
1486 0 : else { return String(id + " must be a vector of integers\n"); }
1487 : }
1488 0 : else{ return String(""); }
1489 : }
1490 :
1491 : // Read Vector<Float> from Record
1492 219 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1493 : {
1494 219 : if( rec.isDefined( id ) )
1495 : {
1496 219 : if( rec.dataType( id )==TpArrayFloat )
1497 : {
1498 73 : rec.get( id , val ); return String("");
1499 : /*
1500 : Array<Float> vec; rec.get(id, vec );
1501 : cout << " vec : " << vec << endl;
1502 : if( vec.shape().nelements()==1 )
1503 : {
1504 : val.resize( vec.shape()[0] );
1505 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec(IPosition(1,i));}
1506 : return String("");
1507 : }
1508 : else { return String(id + " must be a 1D vector of floats"); }
1509 : */
1510 : }
1511 146 : else if ( rec.dataType( id ) ==TpArrayDouble )
1512 : {
1513 0 : Vector<Double> vec; rec.get( id, vec );
1514 0 : val.resize(vec.nelements());
1515 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1516 0 : return String("");
1517 0 : }
1518 146 : else if ( rec.dataType( id ) ==TpArrayInt )
1519 : {
1520 0 : Vector<Int> vec; rec.get( id, vec );
1521 0 : val.resize(vec.nelements());
1522 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1523 0 : return String("");
1524 0 : }
1525 146 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1526 : {
1527 146 : Vector<Bool> vec; rec.get( id, vec );
1528 146 : if( vec.shape().product()==0 ){val.resize(0); return String("");}
1529 0 : else{ return String(id + " must be a vector of strings.\n"); }
1530 : // val.resize(0); return String("");
1531 146 : }
1532 0 : else { return String(id + " must be a vector of floats\n"); }
1533 : }
1534 0 : else{ return String(""); }
1535 : }
1536 :
1537 : // Read Vector<String> from Record
1538 73 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1539 : {
1540 73 : if( rec.isDefined( id ) )
1541 : {
1542 73 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1543 73 : { rec.get( id , val ); return String(""); }
1544 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1545 : {
1546 0 : Vector<Bool> vec; rec.get( id, vec );
1547 0 : if( vec.nelements()==0 ){val.resize(0); return String("");}
1548 0 : else{ return String(id + " must be a vector of strings.\n"); }
1549 0 : }
1550 0 : else { return String(id + " must be a vector of strings.\n");
1551 : }
1552 : }
1553 0 : else{ return String(""); }
1554 : }
1555 :
1556 : // Convert String to Quantity
1557 365 : String SynthesisParams::stringToQuantity(String instr, Quantity& qa) const
1558 : {
1559 : //QuantumHolder qh;
1560 : //String error;
1561 : // if( qh.fromString( error, instr ) ) { qa = qh.asQuantity(); return String(""); }
1562 : //else { return String("Error in converting " + instr + " to a Quantity : " + error + " \n"); }
1563 365 : if ( casacore::Quantity::read( qa, instr ) ) { return String(""); }
1564 0 : else { return String("Error in converting " + instr + " to a Quantity \n"); }
1565 : }
1566 :
1567 : // Convert String to MDirection
1568 : // UR : TODO : If J2000 not specified, make it still work.
1569 2 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1570 : {
1571 : try
1572 : {
1573 2 : String tmpRF, tmpRA, tmpDEC;
1574 2 : std::istringstream iss(instr);
1575 2 : iss >> tmpRF >> tmpRA >> tmpDEC;
1576 2 : if( tmpDEC.length() == 0 )// J2000 may not be specified.
1577 : {
1578 0 : tmpDEC = tmpRA;
1579 0 : tmpRA = tmpRF;
1580 0 : tmpRF = String("J2000");
1581 : }
1582 2 : casacore::Quantity tmpQRA;
1583 2 : casacore::Quantity tmpQDEC;
1584 2 : casacore::Quantity::read(tmpQRA, tmpRA);
1585 2 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1586 :
1587 2 : if(tmpQDEC.getFullUnit()==Unit("deg") && tmpDEC.contains(":")){
1588 0 : LogIO os( LogOrigin("SynthesisParams","stringToMDirection",WHERE) );
1589 : os << LogIO::WARN
1590 : << "You provided the Declination/Latitude value \""<< tmpDEC
1591 : << "\" which is understood to be in units of hours.\n"
1592 : << "If you meant degrees, please replace \":\" by \".\"."
1593 0 : << LogIO::POST;
1594 0 : }
1595 :
1596 : MDirection::Types theRF;
1597 2 : Bool status = MDirection::getType(theRF, tmpRF);
1598 2 : if (!status) {
1599 0 : throw AipsError();
1600 : }
1601 2 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1602 2 : return String("");
1603 2 : }
1604 0 : catch(AipsError &x)
1605 : {
1606 0 : return String("Error in converting '" + instr + "' to MDirection. Need format : 'J2000 19:59:28.500 +40.44.01.50'\n");
1607 0 : }
1608 : }
1609 :
1610 : // Read Quantity from a Record string
1611 365 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1612 : {
1613 365 : if( rec.isDefined( id ) )
1614 : {
1615 292 : if( rec.dataType( id )==TpString )
1616 292 : { String valstr; rec.get( id , valstr ); return stringToQuantity(valstr, val); }
1617 0 : else { return String(id + " must be a string in the format : '1.5GHz' or '0.2arcsec'...'\n"); }
1618 : }
1619 73 : else{ return String(""); }
1620 : }
1621 :
1622 : // Read MDirection from a Record string
1623 75 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1624 : {
1625 75 : if( rec.isDefined( id ) )
1626 : {
1627 2 : if( rec.dataType( id )==TpString )
1628 2 : { String valstr; rec.get( id , valstr ); return stringToMDirection(valstr, val); }
1629 0 : else { return String(id + " must be a string in the format : 'J2000 19:59:28.500 +40.44.01.50'\n"); }
1630 : }
1631 73 : else{ return String(""); }
1632 : }
1633 :
1634 : // Convert MDirection to String
1635 0 : String SynthesisParams::MDirectionToString(MDirection val) const
1636 : {
1637 0 : MVDirection mvpc( val.getAngle() );
1638 0 : MVAngle ra = mvpc.get()(0);
1639 0 : MVAngle dec = mvpc.get()(1);
1640 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1641 0 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1642 0 : }
1643 :
1644 : // Convert Quantity to String
1645 0 : String SynthesisParams::QuantityToString(Quantity val) const
1646 : {
1647 0 : std::ostringstream ss;
1648 : //use digits10 to ensure the conersions involve use full decimal digits guranteed to be
1649 : //correct plus extra digits to deal with least significant digits (or replace with
1650 : // max_digits10 when it is available)
1651 0 : ss.precision(std::numeric_limits<double>::digits10+2);
1652 0 : ss << val;
1653 0 : return ss.str();
1654 : // NOTE - 2017.10.04: It was found (CAS-10773) that we cannot use to_string for this as
1655 : // the decimal place is fixed to 6 digits.
1656 : //TT: change to C++11 to_string which handles double value to string conversion
1657 : //return String(std::to_string( val.getValue(val.getUnit()) )) + val.getUnit() ;
1658 0 : }
1659 :
1660 : // Convert Record contains Quantity or Measure quantities to String
1661 0 : String SynthesisParams::recordQMToString(const Record &rec) const
1662 : {
1663 : Double val;
1664 0 : String unit;
1665 0 : if ( rec.isDefined("m0") )
1666 : {
1667 0 : Record subrec = rec.subRecord("m0");
1668 0 : subrec.get("value",val);
1669 0 : subrec.get("unit",unit);
1670 0 : }
1671 0 : else if (rec.isDefined("value") )
1672 : {
1673 0 : rec.get("value",val);
1674 0 : rec.get("unit",unit);
1675 : }
1676 0 : return String::toString(val) + unit;
1677 0 : }
1678 :
1679 :
1680 : /////////////////////// Selection Parameters
1681 :
1682 219 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1683 : {
1684 219 : setDefaults();
1685 219 : }
1686 :
1687 219 : SynthesisParamsSelect::~SynthesisParamsSelect()
1688 : {
1689 219 : }
1690 :
1691 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1692 0 : operator=(other);
1693 0 : }
1694 :
1695 73 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1696 73 : if(this!=&other) {
1697 73 : msname=other.msname;
1698 73 : spw=other.spw;
1699 73 : freqbeg=other.freqbeg;
1700 73 : freqend=other.freqend;
1701 73 : freqframe=other.freqframe;
1702 73 : field=other.field;
1703 73 : antenna=other.antenna;
1704 73 : timestr=other.timestr;
1705 73 : scan=other.scan;
1706 73 : obs=other.obs;
1707 73 : state=other.state;
1708 73 : uvdist=other.uvdist;
1709 73 : taql=other.taql;
1710 73 : usescratch=other.usescratch;
1711 73 : readonly=other.readonly;
1712 73 : incrmodel=other.incrmodel;
1713 73 : datacolumn=other.datacolumn;
1714 :
1715 : }
1716 73 : return *this;
1717 : }
1718 :
1719 146 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1720 : {
1721 146 : setDefaults();
1722 :
1723 146 : String err("");
1724 :
1725 : try
1726 : {
1727 :
1728 146 : err += readVal( inrec, String("msname"), msname );
1729 :
1730 146 : err += readVal( inrec, String("readonly"), readonly );
1731 146 : err += readVal( inrec, String("usescratch"), usescratch );
1732 :
1733 : // override with entries from savemodel.
1734 146 : String savemodel("");
1735 146 : err += readVal( inrec, String("savemodel"), savemodel );
1736 146 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1737 73 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1738 73 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1739 :
1740 146 : err += readVal( inrec, String("incrmodel"), incrmodel );
1741 :
1742 146 : err += readVal( inrec, String("spw"), spw );
1743 :
1744 : /// -------------------------------------------------------------------------------------
1745 : // Why are these params here ? Repeats of defineimage.
1746 146 : err += readVal( inrec, String("freqbeg"), freqbeg);
1747 146 : err += readVal( inrec, String("freqend"), freqend);
1748 :
1749 146 : String freqframestr( MFrequency::showType(freqframe) );
1750 146 : err += readVal( inrec, String("outframe"), freqframestr);
1751 146 : if( ! MFrequency::getType(freqframe, freqframestr) )
1752 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1753 : /// -------------------------------------------------------------------------------------
1754 :
1755 146 : err += readVal( inrec, String("field"),field);
1756 146 : err += readVal( inrec, String("antenna"),antenna);
1757 146 : err += readVal( inrec, String("timestr"),timestr);
1758 146 : err += readVal( inrec, String("scan"),scan);
1759 146 : err += readVal( inrec, String("obs"),obs);
1760 146 : err += readVal( inrec, String("state"),state);
1761 146 : err += readVal( inrec, String("uvdist"),uvdist);
1762 146 : err += readVal( inrec, String("taql"),taql);
1763 :
1764 146 : err += readVal( inrec, String("datacolumn"),datacolumn);
1765 :
1766 146 : err += verify();
1767 :
1768 146 : }
1769 0 : catch(AipsError &x)
1770 : {
1771 0 : err = err + x.getMesg() + "\n";
1772 0 : }
1773 :
1774 146 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1775 :
1776 146 : }
1777 :
1778 146 : String SynthesisParamsSelect::verify() const
1779 : {
1780 146 : String err;
1781 : // Does the MS exist on disk.
1782 146 : Directory thems( msname );
1783 146 : if( thems.exists() )
1784 : {
1785 : // Is it readable ?
1786 146 : if( ! thems.isReadable() )
1787 0 : { err += "MS " + msname + " is not readable.\n"; }
1788 : // Depending on 'readonly', is the MS writable ?
1789 146 : if( readonly==false && ! thems.isWritable() )
1790 0 : { err += "MS " + msname + " is not writable.\n"; }
1791 : }
1792 : else
1793 0 : { err += "MS does not exist : " + msname + "\n"; }
1794 :
1795 292 : return err;
1796 146 : }
1797 :
1798 365 : void SynthesisParamsSelect::setDefaults()
1799 : {
1800 365 : msname="";
1801 365 : spw="";
1802 365 : freqbeg="";
1803 365 : freqend="";
1804 365 : MFrequency::getType(freqframe,"LSRK");
1805 365 : field="";
1806 365 : antenna="";
1807 365 : timestr="";
1808 365 : scan="";
1809 365 : obs="";
1810 365 : state="";
1811 365 : uvdist="";
1812 365 : taql="";
1813 365 : usescratch=false;
1814 365 : readonly=true;
1815 365 : incrmodel=false;
1816 365 : datacolumn="corrected";
1817 365 : }
1818 :
1819 73 : Record SynthesisParamsSelect::toRecord()const
1820 : {
1821 73 : Record selpar;
1822 73 : selpar.define("msname",msname);
1823 73 : selpar.define("spw",spw);
1824 73 : selpar.define("freqbeg",freqbeg);
1825 73 : selpar.define("freqend",freqend);
1826 73 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1827 : //looks like fromRecord looks for outframe !
1828 73 : selpar.define("outframe", MFrequency::showType(freqframe));
1829 73 : selpar.define("field",field);
1830 73 : selpar.define("antenna",antenna);
1831 73 : selpar.define("timestr",timestr);
1832 73 : selpar.define("scan",scan);
1833 73 : selpar.define("obs",obs);
1834 73 : selpar.define("state",state);
1835 73 : selpar.define("uvdist",uvdist);
1836 73 : selpar.define("taql",taql);
1837 73 : selpar.define("usescratch",usescratch);
1838 73 : selpar.define("readonly",readonly);
1839 73 : selpar.define("incrmodel",incrmodel);
1840 73 : selpar.define("datacolumn",datacolumn);
1841 :
1842 73 : return selpar;
1843 0 : }
1844 :
1845 :
1846 : /////////////////////// Image Parameters
1847 :
1848 219 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1849 : {
1850 219 : setDefaults();
1851 219 : }
1852 :
1853 219 : SynthesisParamsImage::~SynthesisParamsImage()
1854 : {
1855 219 : }
1856 :
1857 146 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1858 146 : if(this != &other){
1859 146 : imageName=other.imageName;
1860 146 : stokes=other.stokes;
1861 146 : startModel.resize(); startModel=other.startModel;
1862 146 : imsize.resize(); imsize=other.imsize;
1863 146 : cellsize.resize(); cellsize=other.cellsize;
1864 146 : projection=other.projection;
1865 146 : useNCP=other.useNCP;
1866 146 : phaseCenter=other.phaseCenter;
1867 146 : phaseCenterFieldId=other.phaseCenterFieldId;
1868 146 : obslocation=other.obslocation;
1869 146 : pseudoi=other.pseudoi;
1870 146 : nchan=other.nchan;
1871 146 : nTaylorTerms=other.nTaylorTerms;
1872 146 : chanStart=other.chanStart;
1873 146 : chanStep=other.chanStep;
1874 146 : freqStart=other.freqStart;
1875 146 : freqStep=other.freqStep;
1876 146 : refFreq=other.refFreq;
1877 146 : velStart=other.velStart;
1878 146 : velStep=other.velStep;
1879 146 : freqFrame=other.freqFrame;
1880 146 : mFreqStart=other.mFreqStart;
1881 146 : mFreqStep=other.mFreqStep;
1882 146 : mVelStart=other.mVelStart;
1883 146 : mVelStep=other.mVelStep;
1884 146 : restFreq.resize(); restFreq=other.restFreq;
1885 146 : start=other.start;
1886 146 : step=other.step;
1887 146 : frame=other.frame;
1888 146 : veltype=other.veltype;
1889 146 : mode=other.mode;
1890 146 : reffreq=other.reffreq;
1891 146 : sysvel=other.sysvel;
1892 146 : sysvelframe=other.sysvelframe;
1893 146 : sysvelvalue=other.sysvelvalue;
1894 146 : qmframe=other.qmframe;
1895 146 : mveltype=other.mveltype;
1896 146 : tststr=other.tststr;
1897 146 : startRecord=other.startRecord;
1898 146 : stepRecord=other.stepRecord;
1899 146 : reffreqRecord=other.reffreqRecord;
1900 146 : sysvelRecord=other.sysvelRecord;
1901 146 : restfreqRecord=other.restfreqRecord;
1902 146 : csysRecord=other.csysRecord;
1903 146 : csys=other.csys;
1904 146 : imshape.resize(); imshape=other.imshape;
1905 146 : freqFrameValid=other.freqFrameValid;
1906 146 : overwrite=other.overwrite;
1907 146 : deconvolver=other.deconvolver;
1908 146 : distance=other.distance;
1909 146 : trackDir=other.trackDir;
1910 146 : trackSource=other.trackSource;
1911 146 : movingSource=other.movingSource;
1912 :
1913 :
1914 :
1915 : }
1916 :
1917 146 : return *this;
1918 :
1919 : }
1920 :
1921 73 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec) {
1922 73 : fromRecord(inrec, False);
1923 73 : }
1924 :
1925 73 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec,
1926 : const casacore::Bool isSingleDish) {
1927 73 : setDefaults();
1928 73 : String err("");
1929 146 : LogIO os(LogOrigin("SynthesisParamsImage", "fromRecord"));
1930 :
1931 : try {
1932 73 : err += readVal(inrec, String("imagename"), imageName);
1933 :
1934 : //// imsize
1935 73 : if (inrec.isDefined("imsize")) {
1936 73 : DataType tp = inrec.dataType("imsize");
1937 73 : if (tp == TpInt) { // A single integer for both dimensions.
1938 : Int npix;
1939 73 : inrec.get("imsize", npix);
1940 73 : imsize.resize(2);
1941 73 : imsize.set(npix);
1942 0 : } else if (tp == TpArrayInt) { // An integer array : [ nx ] or [ nx, ny ]
1943 0 : Vector<Int> ims;
1944 0 : inrec.get("imsize", ims);
1945 0 : if (ims.nelements() == 1) { // [ nx ]
1946 0 : imsize.set(ims[0]);
1947 0 : } else if (ims.nelements() == 2) { // [ nx, ny ]
1948 0 : imsize[0] = ims[0];
1949 0 : imsize[1] = ims[1];
1950 : } else { // [ nx, ny ]
1951 0 : err += "imsize must be either a single integer, or a vector of two integers\n";
1952 : }
1953 0 : } else { // Wrong data type
1954 0 : err += "imsize must be either a single integer, or a vector of two integers\n";
1955 : }
1956 : } //imsize
1957 :
1958 : //// cellsize
1959 73 : if (inrec.isDefined("cell")) {
1960 : try {
1961 73 : DataType tp = inrec.dataType("cell");
1962 73 : if (tp == TpInt ||
1963 73 : tp == TpFloat ||
1964 : tp == TpDouble) {
1965 0 : Double cell = inrec.asDouble("cell");
1966 0 : cellsize.set(Quantity(cell, "arcsec"));
1967 73 : } else if (tp == TpArrayInt ||
1968 73 : tp == TpArrayFloat ||
1969 : tp == TpArrayDouble) {
1970 0 : Vector<Double> cells;
1971 0 : inrec.get("cell", cells);
1972 0 : if (cells.nelements() == 1) { // [ cellx ]
1973 0 : cellsize.set(Quantity(cells[0], "arcsec"));
1974 0 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1975 0 : cellsize[0] = Quantity(cells[0], "arcsec");
1976 0 : cellsize[1] = Quantity(cells[1], "arcsec");
1977 : } else { // Wrong array length
1978 0 : err += "cellsize must be a single integer/string, or a vector of two integers/strings\n";
1979 : }
1980 73 : } else if (tp == TpString) {
1981 73 : String cell;
1982 73 : inrec.get("cell", cell);
1983 73 : Quantity qcell;
1984 73 : err += stringToQuantity(cell, qcell);
1985 73 : cellsize.set(qcell);
1986 73 : } else if (tp == TpArrayString) {
1987 0 : Array<String> cells;
1988 0 : inrec.get("cell", cells);
1989 0 : Vector<String> vcells(cells);
1990 0 : if (cells.nelements() == 1) { // [ cellx ]
1991 0 : Quantity qcell;
1992 0 : err += stringToQuantity(vcells[0], qcell);
1993 0 : cellsize.set(qcell);
1994 0 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1995 0 : err += stringToQuantity(vcells[0], cellsize[0]);
1996 0 : err += stringToQuantity(vcells[1], cellsize[1]);
1997 : } else { // Wrong array length
1998 0 : err += "cellsize must be a single integer/string, or a vector of two integers/strings\n";
1999 : }
2000 0 : } else { // Wrong data type
2001 0 : err += "cellsize must be a single integer/string, or a vector of two integers/strings\n";
2002 : }
2003 0 : } catch (AipsError &x) {
2004 0 : err += "Error reading cellsize : " + x.getMesg();
2005 0 : }
2006 : } // cellsize
2007 :
2008 : //// stokes
2009 73 : err += readVal(inrec, String("stokes"), stokes);
2010 73 : if (stokes.matches("pseudoI")) {
2011 0 : stokes = "I";
2012 0 : pseudoi = true;
2013 : } else {
2014 73 : pseudoi = false;
2015 : }
2016 :
2017 : /// PseudoI
2018 :
2019 : ////nchan
2020 73 : err += readVal(inrec, String("nchan"), nchan);
2021 :
2022 : /// phaseCenter (as a string) . // Add INT support later.
2023 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2024 73 : if (inrec.isDefined("phasecenter")) {
2025 73 : String pcent("");
2026 73 : if (inrec.dataType("phasecenter") == TpString) {
2027 73 : inrec.get("phasecenter", pcent);
2028 73 : if (pcent.length() > 0) { // if it's zero length, it means 'figure out from first field in MS'.
2029 2 : err += readVal(inrec, String("phasecenter"), phaseCenter);
2030 2 : phaseCenterFieldId = -1;
2031 : //// Phase Center Field ID.... if explicitly specified, and not via phasecenter.
2032 : // Need this, to deal with a null phase center being translated to a string to go back out.
2033 2 : err += readVal(inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2034 : } else {
2035 71 : phaseCenterFieldId = -2; // deal with this later in buildCoordinateSystem to assign the first selected field
2036 : }
2037 : }
2038 73 : if (inrec.dataType("phasecenter") == TpInt ||
2039 219 : inrec.dataType("phasecenter") == TpFloat ||
2040 146 : inrec.dataType("phasecenter") == TpDouble) {
2041 : // This will override the previous setting to 0 if the phaseCenter string is zero length.
2042 0 : err += readVal(inrec, String("phasecenter"), phaseCenterFieldId);
2043 : }
2044 73 : if ((inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter") != TpInt
2045 73 : && inrec.dataType("phasecenter") != TpFloat && inrec.dataType("phasecenter") != TpDouble)) {
2046 0 : err += String("Cannot set phasecenter. Please specify a string or int\n");
2047 : }
2048 73 : }
2049 :
2050 : //// Projection
2051 73 : if (inrec.isDefined("projection")) {
2052 73 : if (inrec.dataType("projection") == TpString) {
2053 73 : String pstr;
2054 73 : inrec.get("projection", pstr);
2055 : try {
2056 73 : if (pstr.matches("NCP")) {
2057 0 : pstr = "SIN";
2058 0 : useNCP = true;
2059 : }
2060 73 : projection = Projection::type(pstr);
2061 0 : } catch (AipsError &x) {
2062 0 : err += String("Invalid projection code : " + pstr);
2063 0 : }
2064 73 : } else {
2065 0 : err += "projection must be a string\n";
2066 : }
2067 : } //projection
2068 :
2069 : // Frequency frame stuff.
2070 73 : err += readVal(inrec, String("specmode"), mode);
2071 : // Alias for 'mfs' is 'cont'
2072 73 : if (mode == "cont") mode = "mfs";
2073 :
2074 73 : err += readVal(inrec, String("outframe"), frame);
2075 73 : qmframe = "";
2076 : // mveltype is only set when start/step is given in mdoppler
2077 73 : mveltype = "";
2078 : //start
2079 73 : String startType("");
2080 73 : String widthType("");
2081 73 : if (inrec.isDefined("start")) {
2082 73 : if (inrec.dataType("start") == TpInt) {
2083 0 : err += readVal(inrec, String("start"), chanStart);
2084 0 : start = String::toString(chanStart);
2085 0 : startType = "chan";
2086 73 : } else if (inrec.dataType("start") == TpString) {
2087 73 : err += readVal(inrec, String("start"), start);
2088 73 : if (start.contains("Hz")) {
2089 0 : stringToQuantity(start, freqStart);
2090 0 : startType = "freq";
2091 73 : } else if (start.contains("m/s")) {
2092 0 : stringToQuantity(start, velStart);
2093 0 : startType = "vel";
2094 : }
2095 0 : } else if (inrec.dataType("start") == TpRecord) {
2096 : //record can be freq in Quantity or MFreaquency or vel in Quantity or
2097 : //MRadialVelocity or Doppler (by me.todoppler())
2098 : // ** doppler => converted to radialvel with frame
2099 0 : startRecord = inrec.subRecord("start");
2100 0 : if (startRecord.isDefined("m0")) {
2101 : //must be a measure
2102 0 : String mtype;
2103 0 : startRecord.get("type", mtype);
2104 0 : if (mtype == "frequency") {
2105 : //mfrequency
2106 0 : startRecord.get(String("refer"), qmframe);
2107 0 : if (frame != "" && frame != qmframe) {
2108 : // should emit warning to the logger
2109 0 : cerr << "The frame in start:" << qmframe << " Override frame=" << frame << endl;
2110 : }
2111 0 : start = recordQMToString(startRecord);
2112 0 : stringToQuantity(start, freqStart);
2113 0 : startType = "freq";
2114 0 : } else if (mtype == "radialvelocity") {
2115 : //mradialvelocity
2116 0 : startRecord.get(String("refer"), qmframe);
2117 0 : if (frame != "" && frame != qmframe) {
2118 : // should emit warning to the logger
2119 0 : cerr << "The frame in start:" << qmframe << " Override frame=" << frame << endl;
2120 : }
2121 0 : start = recordQMToString(startRecord);
2122 0 : stringToQuantity(start, velStart);
2123 0 : startType = "vel";
2124 0 : } else if (mtype == "doppler") {
2125 : //use veltype in mdoppler
2126 : //start = MDopToVelString(startRecord);
2127 0 : start = recordQMToString(startRecord);
2128 0 : stringToQuantity(start, velStart);
2129 0 : startRecord.get(String("refer"), mveltype);
2130 0 : mveltype.downcase();
2131 0 : startType = "vel";
2132 : }
2133 0 : } else {
2134 0 : start = recordQMToString(startRecord);
2135 0 : if (start.contains("Hz")) {
2136 0 : stringToQuantity(start, freqStart);
2137 0 : startType = "freq";
2138 0 : } else if (start.contains("m/s")) {
2139 0 : stringToQuantity(start, velStart);
2140 0 : startType = "vel";
2141 : } else {
2142 0 : err += String("Unrecognized Quantity unit for start, must contain m/s or Hz\n");
2143 : }
2144 : }
2145 : } else {
2146 0 : err += String("start must be an integer, a string, or frequency/velocity in Quantity/Measure\n");
2147 : }
2148 : }
2149 :
2150 : //step
2151 73 : if (inrec.isDefined("width")) {
2152 73 : if (inrec.dataType("width") == TpInt) {
2153 0 : err += readVal(inrec, String("width"), chanStep);
2154 0 : step = String::toString(chanStep);
2155 0 : widthType = "chan";
2156 73 : } else if (inrec.dataType("width") == TpString) {
2157 73 : err += readVal(inrec, String("width"), step);
2158 73 : if (step.contains("Hz")) {
2159 0 : stringToQuantity(step, freqStep);
2160 0 : widthType = "freq";
2161 73 : } else if (step.contains("m/s")) {
2162 0 : stringToQuantity(step, velStep);
2163 0 : widthType = "vel";
2164 : }
2165 0 : } else if (inrec.dataType("width") == TpRecord) {
2166 : //record can be freq in Quantity or MFreaquency or vel in Quantity or
2167 : //MRadialVelocity or Doppler (by me.todoppler())
2168 : // ** doppler => converted to radialvel with frame
2169 0 : stepRecord = inrec.subRecord("width");
2170 0 : if (stepRecord.isDefined("m0")) {
2171 : //must be a measure
2172 0 : String mtype;
2173 0 : stepRecord.get("type", mtype);
2174 0 : if (mtype == "frequency") {
2175 : //mfrequency
2176 0 : stepRecord.get(String("refer"), qmframe);
2177 0 : if (frame != "" && frame != qmframe) {
2178 : // should emit warning to the logger
2179 0 : cerr << "The frame in step:" << qmframe << " Override frame=" << frame << endl;
2180 : }
2181 0 : step = recordQMToString(stepRecord);
2182 0 : stringToQuantity(step, freqStep);
2183 0 : widthType = "freq";
2184 0 : } else if (mtype == "radialvelocity") {
2185 : //mradialvelocity
2186 0 : stepRecord.get(String("refer"), qmframe);
2187 0 : if (frame != "" && frame != qmframe) {
2188 : // should emit warning to the logger
2189 0 : cerr << "The frame in step:" << qmframe << " Override frame=" << frame << endl;
2190 : }
2191 0 : step = recordQMToString(stepRecord);
2192 0 : stringToQuantity(step, velStep);
2193 0 : widthType = "vel";
2194 0 : } else if (mtype == "doppler") {
2195 : //step = MDopToVelString(stepRecord);
2196 0 : step = recordQMToString(stepRecord);
2197 0 : stringToQuantity(step, velStep);
2198 0 : startRecord.get(String("refer"), mveltype);
2199 0 : mveltype.downcase();
2200 0 : widthType = "vel";
2201 : }
2202 0 : } else {
2203 0 : step = recordQMToString(stepRecord);
2204 0 : if (step.contains("Hz")) {
2205 0 : stringToQuantity(step, freqStep);
2206 0 : widthType = "freq";
2207 0 : } else if (step.contains("m/s")) {
2208 0 : stringToQuantity(step, velStep);
2209 0 : widthType = "vel";
2210 : } else {
2211 0 : err += String("Unrecognized Quantity unit for step, must contain m/s or Hz\n");
2212 : }
2213 : }
2214 : } else {
2215 0 : err += String("step must be an integer, a string, or frequency/velocity in Quantity/Measure\n");
2216 : }
2217 : }
2218 :
2219 : //check for start, width unit consistentcy
2220 73 : if (startType != widthType && startType != "" && widthType != "") {
2221 0 : err += String("Cannot mix start and width with different unit types (e.g. km/s vs. Hz)\n");
2222 : }
2223 :
2224 : //reffreq (String, Quantity, or Measure)
2225 73 : if (inrec.isDefined("reffreq")) {
2226 73 : if (inrec.dataType("reffreq") == TpString) {
2227 73 : err += readVal(inrec, String("reffreq"), refFreq);
2228 0 : } else if (inrec.dataType("reffreq") == TpRecord) {
2229 0 : String reffreqstr;
2230 0 : reffreqRecord = inrec.subRecord("reffreq");
2231 0 : if (reffreqRecord.isDefined("m0")) {
2232 0 : String mtype;
2233 0 : reffreqRecord.get("type", mtype);
2234 0 : if (mtype == "frequency") {
2235 0 : reffreqstr = recordQMToString(reffreqRecord);
2236 0 : stringToQuantity(reffreqstr, refFreq);
2237 : } else {
2238 0 : err += String("Unrecognized Measure for reffreq, must be a frequency measure\n");
2239 : }
2240 0 : } else {
2241 0 : reffreqstr = recordQMToString(reffreqRecord);
2242 0 : if (reffreqstr.contains("Hz")) {
2243 0 : stringToQuantity(reffreqstr, refFreq);
2244 : } else {
2245 0 : err += String("Unrecognized Quantity unit for reffreq, must contain Hz\n");
2246 : }
2247 : }
2248 0 : } else {
2249 0 : err += String("reffreq must be a string, or frequency in Quantity/Measure\n");
2250 : }
2251 : }
2252 :
2253 73 : err += readVal(inrec, String("veltype"), veltype);
2254 73 : veltype = mveltype != "" ? mveltype : veltype;
2255 : // sysvel (String, Quantity)
2256 73 : if (inrec.isDefined("sysvel")) {
2257 73 : if (inrec.dataType("sysvel") == TpString) {
2258 73 : err += readVal(inrec, String("sysvel"), sysvel);
2259 0 : } else if (inrec.dataType("sysvel") == TpRecord) {
2260 0 : sysvelRecord = inrec.subRecord("sysvel");
2261 0 : sysvel = recordQMToString(sysvelRecord);
2262 0 : if (sysvel == "" || !sysvel.contains("m/s")) {
2263 0 : err += String("Unrecognized Quantity unit for sysvel, must contain m/s\n");
2264 : }
2265 : } else {
2266 0 : err += String("sysvel must be a string, or velocity in Quantity\n");
2267 : }
2268 : }
2269 73 : err += readVal(inrec, String("sysvelframe"), sysvelframe);
2270 :
2271 : // rest frequencies (record or vector of Strings)
2272 73 : if (inrec.isDefined("restfreq")) {
2273 73 : Vector<String> rfreqs(0);
2274 73 : Record restfreqSubRecord;
2275 73 : if (inrec.dataType("restfreq") == TpRecord) {
2276 0 : restfreqRecord = inrec.subRecord("restfreq");
2277 : // assume multiple restfreqs are index as '0','1'..
2278 0 : if (restfreqRecord.isDefined("0")) {
2279 0 : rfreqs.resize(restfreqRecord.nfields());
2280 0 : for (uInt fr = 0; fr < restfreqRecord.nfields(); fr++) {
2281 0 : restfreqSubRecord = restfreqRecord.subRecord(String::toString(fr));
2282 0 : rfreqs[fr] = recordQMToString(restfreqSubRecord);
2283 : }
2284 : }
2285 73 : } else if (inrec.dataType("restfreq") == TpArrayString) {
2286 0 : err += readVal(inrec, String("restfreq"), rfreqs);
2287 73 : } else if (inrec.dataType("restfreq") == TpString) {
2288 0 : rfreqs.resize(1);
2289 0 : err += readVal(inrec, String("restfreq"), rfreqs[0]);
2290 : }
2291 73 : restFreq.resize(rfreqs.nelements());
2292 73 : for (uInt fr = 0; fr < rfreqs.nelements(); fr++) {
2293 0 : err += stringToQuantity(rfreqs[fr], restFreq[fr]);
2294 : }
2295 73 : } // if def restfreq
2296 :
2297 : // optional - coordsys, imshape
2298 : // if exist use them. May need a consistency check with the rest of impars?
2299 73 : if (inrec.isDefined("csys")) {
2300 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
2301 0 : if (inrec.dataType("csys") == TpRecord) {
2302 : //csysRecord = inrec.subRecord("csys");
2303 0 : csysRecord.defineRecord("coordsys", inrec.subRecord("csys"));
2304 : }
2305 0 : if (inrec.isDefined("imshape")) {
2306 0 : if (inrec.dataType("imshape") == TpArrayInt) {
2307 0 : err += readVal(inrec, String("imshape"), imshape);
2308 : }
2309 : }
2310 : }
2311 :
2312 : //String freqframestr( MFrequency::showType(freqFrame) );
2313 : //err += readVal( inrec, String("outframe"), freqframestr);
2314 : //if( ! MFrequency::getType(freqFrame, freqframestr) )
2315 : // { err += "Invalid Frequency Frame " + freqframestr ; }
2316 :
2317 73 : String freqframestr = (frame != "" && qmframe != "") ? qmframe : frame;
2318 73 : if (frame != "" && !MFrequency::getType(freqFrame, freqframestr)) {
2319 0 : err += "Invalid Frequency Frame " + freqframestr;
2320 : }
2321 73 : err += readVal(inrec, String("restart"), overwrite);
2322 :
2323 73 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2324 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2325 73 : if (inrec.isDefined("startmodel")) {
2326 73 : if (inrec.dataType("startmodel") == TpString) {
2327 73 : String onemodel;
2328 73 : err += readVal(inrec, String("startmodel"), onemodel);
2329 73 : if (onemodel.length() > 0) {
2330 0 : startModel.resize(1);
2331 0 : startModel[0] = onemodel;
2332 : } else {
2333 73 : startModel.resize();
2334 : }
2335 73 : } else if (inrec.dataType("startmodel") == TpArrayString ||
2336 0 : inrec.dataType("startmodel") == TpArrayBool) {
2337 0 : err += readVal(inrec, String("startmodel"), startModel);
2338 : } else {
2339 0 : err += String("startmodel must be either a string(singleterm) or a list of strings(multiterm)\n");
2340 : }
2341 : }
2342 :
2343 73 : err += readVal(inrec, String("nterms"), nTaylorTerms);
2344 73 : err += readVal(inrec, String("deconvolver"), deconvolver);
2345 :
2346 : // Force nchan=1 for anything other than cube modes...
2347 73 : if (mode == "mfs") nchan = 1;
2348 : //read obslocation
2349 73 : if (inrec.isDefined("obslocation_rec")) {
2350 0 : String errorobs;
2351 0 : const Record obsrec = inrec.asRecord("obslocation_rec");
2352 0 : MeasureHolder mh;
2353 0 : if (!mh.fromRecord(errorobs, obsrec)) {
2354 0 : err += errorobs;
2355 : }
2356 0 : obslocation = mh.asMPosition();
2357 0 : }
2358 73 : err += verify(isSingleDish);
2359 73 : }
2360 0 : catch (AipsError &x) {
2361 0 : err = err + x.getMesg() + "\n";
2362 0 : }
2363 :
2364 73 : if (err.length() > 0) throw(AipsError("Invalid Image Parameter set : " + err));
2365 73 : }
2366 :
2367 0 : String SynthesisParamsImage::MDopToVelString(Record &rec)
2368 : {
2369 0 : if( rec.isDefined("type") ) {
2370 0 : String measType;
2371 0 : String unit;
2372 0 : Double val = 0;
2373 0 : rec.get("type", measType);
2374 0 : if(measType=="doppler") {
2375 0 : rec.get(String("refer"), mveltype);
2376 0 : Record dopRecord = rec.subRecord("m0");
2377 0 : String dopstr = recordQMToString(dopRecord);
2378 : //cerr<<"dopstr="<<dopstr<<endl;
2379 : MRadialVelocity::Types mvType;
2380 : //use input frame
2381 0 : qmframe = frame!=""? frame: "LSRK";
2382 0 : MRadialVelocity::getType(mvType, qmframe);
2383 : MDoppler::Types mdType;
2384 0 : MDoppler::getType(mdType, mveltype);
2385 0 : MDoppler dop(Quantity(val,unit), mdType);
2386 0 : MRadialVelocity mRadVel(MRadialVelocity::fromDoppler(dop, mvType));
2387 0 : Double velval = mRadVel.get("m/s").getValue();
2388 0 : return start = String::toString(velval) + String("m/s");
2389 0 : }
2390 0 : else { return String("");}
2391 0 : }
2392 0 : else { return String("");}
2393 : }
2394 :
2395 0 : String SynthesisParamsImage::verify() const
2396 : {
2397 0 : verify(False);
2398 0 : }
2399 :
2400 73 : String SynthesisParamsImage::verify(const casacore::Bool isSingleDish) const
2401 : {
2402 73 : LogIO os;
2403 73 : String err;
2404 :
2405 73 : if(imageName=="") {err += "Please supply an image name\n";}
2406 :
2407 73 : if(imsize.nelements() != 2){err += "imsize must be a vector of 2 Ints\n";}
2408 73 : if(cellsize.nelements() != 2) {err += "cellsize must be a vector of 2 Quantities\n";}
2409 73 : if(cellsize[0].getValue() == 0.0 || cellsize[1].getValue() == 0.0){
2410 0 : err += "cellsize must be nonzero\n";
2411 : }
2412 :
2413 : //// default is nt=2 but deconvolver != mtmfs by default.
2414 : // if( nchan>1 and nTaylorTerms>1 )
2415 : // {err += "Cannot have more than one channel with ntaylorterms>1\n";}
2416 :
2417 73 : if((mode=="mfs") && nchan > 1){
2418 0 : err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";
2419 : }
2420 :
2421 73 : if(! stokes.matches("I") && ! stokes.matches("Q") &&
2422 0 : ! stokes.matches("U") && ! stokes.matches("V") &&
2423 0 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2424 0 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2425 0 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2426 0 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2427 0 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2428 73 : ! stokes.matches("IQU") && ! stokes.matches("IUV") &&
2429 0 : ! stokes.matches("IQUV")){
2430 :
2431 0 : err += "Stokes " + stokes + " is an unsupported option \n";
2432 : }
2433 :
2434 : /// err += verifySpectralSetup();
2435 :
2436 : // Allow only one starting model. No additions to be done.
2437 73 : if(startModel.nelements()>0){
2438 0 : if(deconvolver!="mtmfs"){
2439 :
2440 0 : if( startModel.nelements()!=1 ){
2441 0 : err += String("Only one startmodel image is allowed.\n");
2442 : }
2443 : else{
2444 0 : File fp(imageName+String(".model"));
2445 0 : if(fp.exists()) err += "Model " + imageName+".model exists, but a starting model of " + startModel[0] + " is also being requested. Please either reset startmodel='' to use what already exists, or delete " + imageName + ".model so that it uses the new model specified in startmodel.";
2446 0 : }
2447 : }
2448 : else{// mtmfs
2449 0 : File fp(imageName+String(".model.tt0"));
2450 0 : if(fp.exists()){
2451 0 : err += "Model " + imageName+".model.tt* exists, but a starting model of ";
2452 0 : for (uInt i=0;i<startModel.nelements();i++){err += startModel[i] + ",";}
2453 0 : err +=" is also being requested. Please either reset startmodel='' to use what already exists, or delete " + imageName + ".model.tt* so that it uses the new model specified in startmodel";
2454 : }
2455 0 : }
2456 :
2457 : // Check that startmodel exists on disk !
2458 0 : for(uInt ss=0;ss<startModel.nelements();ss++){
2459 0 : File fp( startModel[ss] );
2460 0 : if(!fp.exists()){
2461 0 : err += "Startmodel " + startModel[ss] + " cannot be found on disk.";
2462 : }
2463 0 : }
2464 : }
2465 :
2466 : /// Check imsize for efficiency.
2467 73 : Int imxnew = SynthesisUtilMethods::getOptimumSize(imsize[0]);
2468 73 : Int imynew = SynthesisUtilMethods::getOptimumSize(imsize[1]);
2469 73 : if((imxnew != imsize[0] || imynew != imsize[1]) && isSingleDish == False){
2470 0 : LogIO os(LogOrigin("SynthesisParamsImage","fromRecord",WHERE));
2471 0 : if( imxnew != imsize[0] ) {
2472 0 : os << LogIO::WARN << "imsize with "+String::toString(imsize[0])+
2473 0 : " pixels is not an efficient imagesize. Try "+
2474 0 : String::toString(imxnew)+" instead." << LogIO::POST;
2475 : }
2476 0 : if( imsize[0] != imsize[1] && imynew != imsize[1] ){
2477 0 : os << LogIO::WARN << "imsize with "+String::toString(imsize[1])+
2478 0 : " pixels is not an efficient imagesize. Try "+
2479 0 : String::toString(imynew)+" instead." << LogIO::POST;
2480 : }
2481 0 : }
2482 146 : return err;
2483 73 : }// verify()
2484 :
2485 : /*
2486 : // Convert all user options to LSRK freqStart, freqStep,
2487 : // Could have (optional) log messages coming out of this function, to tell the user what the
2488 : // final frequency setup is ?
2489 :
2490 : String SynthesisParamsImage::verifySpectralSetup()
2491 : {
2492 : }
2493 : */
2494 :
2495 292 : void SynthesisParamsImage::setDefaults()
2496 : {
2497 : // Image definition parameters
2498 292 : imageName = String("");
2499 292 : imsize.resize(2); imsize.set(100);
2500 292 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2501 292 : stokes="I";
2502 292 : phaseCenter=MDirection();
2503 292 : phaseCenterFieldId=-1;
2504 292 : projection=Projection::SIN;
2505 292 : useNCP=false;
2506 292 : startModel=Vector<String>(0);
2507 292 : freqFrameValid=True;
2508 292 : overwrite=false;
2509 : // PseudoI
2510 292 : pseudoi=false;
2511 :
2512 : // Spectral coordinates
2513 292 : nchan=1;
2514 292 : mode="mfs";
2515 292 : start="";
2516 292 : step="";
2517 292 : chanStart=0;
2518 292 : chanStep=1;
2519 : //freqStart=Quantity(0,"Hz");
2520 : //freqStep=Quantity(0,"Hz");
2521 : //velStart=Quantity(0,"m/s");
2522 : //velStep=Quantity(0,"m/s");
2523 292 : freqStart=Quantity(0,"");
2524 292 : freqStep=Quantity(0,"");
2525 292 : velStart=Quantity(0,"");
2526 292 : velStep=Quantity(0,"");
2527 292 : veltype=String("radio");
2528 292 : restFreq.resize(0);
2529 292 : refFreq = Quantity(0,"Hz");
2530 292 : frame = "";
2531 292 : freqFrame=MFrequency::LSRK;
2532 292 : sysvel="";
2533 292 : sysvelframe="";
2534 292 : sysvelvalue=Quantity(0.0,"m/s");
2535 292 : nTaylorTerms=1;
2536 292 : deconvolver="hogbom";
2537 : ///csysRecord=Record();
2538 292 : }
2539 :
2540 0 : Record SynthesisParamsImage::toRecord() const
2541 : {
2542 0 : Record impar;
2543 0 : impar.define("imagename", imageName);
2544 0 : impar.define("imsize", imsize);
2545 0 : Vector<String> cells(2);
2546 0 : cells[0] = QuantityToString( cellsize[0] );
2547 0 : cells[1] = QuantityToString( cellsize[1] );
2548 0 : impar.define("cell", cells );
2549 0 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2550 0 : else{impar.define("stokes", stokes);}
2551 0 : impar.define("nchan", nchan);
2552 0 : impar.define("nterms", nTaylorTerms);
2553 0 : impar.define("deconvolver",deconvolver);
2554 0 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2555 0 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2556 0 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2557 :
2558 0 : impar.define("specmode", mode );
2559 : // start and step can be one of these types
2560 0 : if( start!="" )
2561 : {
2562 0 : if( !start.contains("Hz") && !start.contains("m/s") &&
2563 0 : String::toInt(start) == chanStart )
2564 : {
2565 0 : impar.define("start",chanStart);
2566 : }
2567 0 : else if( startRecord.nfields() > 0 )
2568 : {
2569 0 : impar.defineRecord("start", startRecord );
2570 : }
2571 : else
2572 : {
2573 0 : impar.define("start",start);
2574 : }
2575 : }
2576 : else {
2577 0 : impar.define("start", start );
2578 : }
2579 0 : if( step!="" )
2580 : {
2581 0 : if( !step.contains("Hz") && !step.contains("m/s") &&
2582 0 : String::toInt(step) == chanStep )
2583 : {
2584 0 : impar.define("width", chanStep);
2585 : }
2586 0 : else if( stepRecord.nfields() > 0 )
2587 : {
2588 0 : impar.defineRecord("width",stepRecord);
2589 : }
2590 : else
2591 : {
2592 0 : impar.define("width",step);
2593 : }
2594 : }
2595 : else
2596 : {
2597 0 : impar.define("width", step);
2598 : }
2599 : //impar.define("chanstart", chanStart );
2600 : //impar.define("chanstep", chanStep );
2601 : //impar.define("freqstart", QuantityToString( freqStart ));
2602 : //impar.define("freqstep", QuantityToString( freqStep ) );
2603 : //impar.define("velstart", QuantityToString( velStart ));
2604 : //impar.define("velstep", QuantityToString( velStep ) );
2605 0 : impar.define("veltype", veltype);
2606 0 : if (restfreqRecord.nfields() != 0 )
2607 : {
2608 0 : impar.defineRecord("restfreq", restfreqRecord);
2609 : }
2610 : else
2611 : {
2612 0 : Vector<String> rfs( restFreq.nelements() );
2613 0 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2614 0 : impar.define("restfreq", rfs);
2615 0 : }
2616 : //impar.define("reffreq", QuantityToString(refFreq));
2617 : //reffreq
2618 0 : if( reffreqRecord.nfields() != 0 )
2619 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2620 : else
2621 0 : { impar.define("reffreq",reffreq); }
2622 : //impar.define("reffreq", reffreq );
2623 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2624 0 : impar.define("outframe", frame );
2625 : //sysvel
2626 0 : if( sysvelRecord.nfields() != 0 )
2627 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2628 : else
2629 0 : { impar.define("sysvel", sysvel );}
2630 0 : impar.define("sysvelframe", sysvelframe );
2631 :
2632 0 : impar.define("restart",overwrite );
2633 0 : impar.define("freqframevalid", freqFrameValid);
2634 0 : impar.define("startmodel", startModel );
2635 :
2636 0 : if( csysRecord.isDefined("coordsys") )
2637 : {
2638 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2639 0 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2640 0 : impar.define("imshape", imshape);
2641 : }
2642 : // else cout << " NO CSYS INFO to write to output record " << endl;
2643 : ///Now save obslocation
2644 0 : Record tmprec;
2645 0 : String err;
2646 0 : MeasureHolder mh(obslocation);
2647 0 : if(mh.toRecord(err, tmprec)){
2648 0 : impar.defineRecord("obslocation_rec", tmprec);
2649 : }
2650 : else{
2651 0 : throw(AipsError("failed to save obslocation to record"));
2652 : }
2653 0 : return impar;
2654 0 : }
2655 :
2656 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2657 : /////////////////////////// Build a coordinate system. ////////////////////////////////////////
2658 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2659 : //// To be used from SynthesisImager, to construct the images it needs
2660 : //// To also be connected to a 'makeimage' method of the synthesisimager tool.
2661 : //// ( need to supply MS only to add 'ObsInfo' to the csys )
2662 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2663 :
2664 :
2665 :
2666 73 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystem(vi::VisibilityIterator2& vi2, const std::map<Int, std::map<Int, Vector<Int> > >& chansel, Block<const MeasurementSet *> mss)
2667 :
2668 : {
2669 : //vi2.getImpl()->spectralWindows( spwids );
2670 : //The above is not right
2671 : //////////// ///Kludge to find all spw selected
2672 : //std::vector<Int> pushspw;
2673 73 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2674 73 : vi2.originChunks();
2675 73 : vi2.origin();
2676 : /// This version uses the new vi2/vb2
2677 : // get the first ms for multiple MSes
2678 : //MeasurementSet msobj=vi2.ms();
2679 73 : Int fld=vb->fieldId()(0);
2680 :
2681 : //handling first ms only
2682 73 : Double gfreqmax=-1.0;
2683 73 : Double gdatafend=-1.0;
2684 73 : Double gdatafstart=1e14;
2685 73 : Double gfreqmin=1e14;
2686 73 : Vector<Int> spwids0;
2687 73 : Int j=0;
2688 73 : Int minfmsid=0;
2689 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2690 73 : Double imStartFreq=getCubeImageStartFreq();
2691 73 : std::vector<Int> sourceMsWithStartFreq;
2692 :
2693 :
2694 146 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2695 : //auto forMS0=chansel.find(0);
2696 73 : map<Int, Vector<Int> > spwsels=forMS0->second;
2697 73 : Int nspws=spwsels.size();
2698 73 : Vector<Int> spwids(nspws);
2699 73 : Vector<Int> nChannels(nspws);
2700 73 : Vector<Int> firstChannels(nspws);
2701 : //Vector<Int> channelIncrement(nspws);
2702 :
2703 73 : Int k=0;
2704 146 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2705 73 : spwids[k]=it->first;
2706 73 : nChannels[k]=(it->second)[0];
2707 73 : firstChannels[k]=(it->second)[1];
2708 : }
2709 73 : if(j==0) {
2710 73 : spwids0.resize();
2711 73 : spwids0=spwids;
2712 : }
2713 : // std::tie (spwids, nChannels, firstChannels, channelIncrement)=(static_cast<vi::VisibilityIteratorImpl2 * >(vi2.getImpl()))->getChannelInformation(false);
2714 :
2715 : //cerr << "SPWIDS "<< spwids << " nchan " << nChannels << " firstchan " << firstChannels << endl;
2716 :
2717 : //////////////////This returns junk for multiple ms CAS-9994..so kludged up along with spw kludge
2718 : //Vector<Int> flds;
2719 : //vi2.getImpl()->fieldIds( flds );
2720 : //AlwaysAssert( flds.nelements()>0 , AipsError );
2721 : //fld = flds[0];
2722 73 : Double freqmin=0, freqmax=0;
2723 73 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2724 :
2725 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2726 73 : MFrequency::Types dataFrame=(MFrequency::Types)MSColumns(*mss[j]).spectralWindow().measFreqRef()(spwids[0]);
2727 :
2728 : Double datafstart, datafend;
2729 : //VisBufferUtil::getFreqRange(datafstart, datafend, vi2, dataFrame );
2730 : //cerr << std::setprecision(12) << "before " << datafstart << " " << datafend << endl;
2731 73 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2732 : //cerr << "after " << datafstart << " " << datafend << endl;
2733 73 : if((datafstart > datafend) || !status)
2734 0 : throw(AipsError("spw selection failed"));
2735 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2736 :
2737 73 : if (mode=="cubedata") {
2738 0 : freqmin = datafstart;
2739 0 : freqmax = datafend;
2740 : }
2741 73 : else if(mode == "cubesource"){
2742 0 : if(!trackSource){
2743 0 : throw(AipsError("Cannot be in cubesource without tracking a moving source"));
2744 : }
2745 0 : String ephemtab(movingSource);
2746 0 : if(movingSource=="TRACKFIELD"){
2747 0 : Int fieldID=MSColumns(*mss[j]).fieldId()(0);
2748 0 : ephemtab=Path(MSColumns(*mss[j]).field().ephemPath(fieldID)).absoluteName();
2749 : }
2750 0 : MEpoch refep=MSColumns(*mss[j]).timeMeas()(0);
2751 0 : Quantity refsysvel;
2752 0 : MSUtil::getFreqRangeAndRefFreqShift(freqmin,freqmax,refsysvel, refep, spwids,firstChannels, nChannels, *mss[j], ephemtab, trackDir, true);
2753 0 : if(j==0)
2754 0 : sysvelvalue=refsysvel;
2755 : /*Double freqMinTopo, freqMaxTopo;
2756 : MSUtil::getFreqRangeInSpw( freqMinTopo, freqMaxTopo, spwids, firstChannels,
2757 : nChannels,*mss[j], freqFrameValid? MFrequency::TOPO:MFrequency::REST , True);
2758 : cerr << std::setprecision(10) << (freqmin-freqMinTopo) << " " << (freqmax-freqMaxTopo) << endl;
2759 : sysfreqshift=((freqmin-freqMinTopo)+(freqmax-freqMaxTopo))/2.0;
2760 : */
2761 0 : }
2762 : else {
2763 : //VisBufferUtil::getFreqRange(freqmin,freqmax, vi2, freqFrameValid? freqFrame:MFrequency::REST );
2764 : //cerr << "before " << freqmin << " " << freqmax << endl;
2765 146 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2766 146 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2767 : //cerr << "after " << freqmin << " " << freqmax << endl;
2768 : }
2769 :
2770 :
2771 :
2772 :
2773 73 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2774 73 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2775 73 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2776 73 : if(datafend > gdatafend) gdatafend=datafend;
2777 : // pick ms to use for setting spectral coord for output images
2778 : // when startfreq is specified find first ms that it fall within the freq range
2779 : // of the ms (with channel selection applied).
2780 : // startfreq is converted to the data frame freq based on Measure ref (for the direction, epech, location)
2781 : // of that ms.
2782 73 : if(imStartFreq > 0.0 && imStartFreq >= freqmin && imStartFreq <= freqmax){
2783 0 : if(mode != "cubesource"){
2784 0 : minfmsid=j;
2785 0 : spwids0.resize();
2786 0 : spwids0=spwids;
2787 0 : vi2.originChunks();
2788 0 : vi2.origin();
2789 0 : while(vb->msId() != j && vi2.moreChunks() ){
2790 0 : vi2.nextChunk();
2791 0 : vi2.origin();
2792 : }
2793 0 : fld=vb->fieldId()(0);
2794 :
2795 : }
2796 : else{
2797 0 : sourceMsWithStartFreq.push_back(j);
2798 : }
2799 : }
2800 :
2801 73 : }
2802 73 : if(sourceMsWithStartFreq.size() > 1){
2803 0 : auto result = std::find(std::begin(sourceMsWithStartFreq), std::end(sourceMsWithStartFreq), 0);
2804 0 : if(result == std::end(sourceMsWithStartFreq)){
2805 0 : throw(AipsError("Reorder the input list of MSs so that MS "+String::toString( sourceMsWithStartFreq[0])+ "is first to match startfreq you provided"));
2806 : }
2807 : }
2808 73 : MeasurementSet msobj = *mss[minfmsid];
2809 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2810 146 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2811 73 : }
2812 :
2813 :
2814 0 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystem(ROVisibilityIterator* rvi )
2815 : {
2816 : /// This version uses the old vi/vb
2817 : // get the first ms for multiple MSes
2818 0 : MeasurementSet msobj=rvi->getMeasurementSet();
2819 0 : Vector<Int> spwids;
2820 0 : Vector<Int> nvischan;
2821 0 : rvi->allSelectedSpectralWindows(spwids,nvischan);
2822 0 : Int fld = rvi->fieldId();
2823 0 : Double freqmin=0, freqmax=0;
2824 : Double datafstart, datafend;
2825 : //freqFrameValid=(freqFrame != MFrequency::REST || mode != "cubedata" );
2826 0 : freqFrameValid=(freqFrame != MFrequency::REST );
2827 0 : MSColumns msc(msobj);
2828 0 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2829 0 : rvi->getFreqInSpwRange(datafstart, datafend, dataFrame );
2830 0 : if (mode=="cubedata") {
2831 0 : freqmin = datafstart;
2832 0 : freqmax = datafend;
2833 : }
2834 : else {
2835 0 : rvi->getFreqInSpwRange(freqmin,freqmax,freqFrameValid? freqFrame:MFrequency::REST );
2836 : }
2837 : // Following three lines are kind of redundant but need to get freq range in the data frame to be used
2838 : // to select channel range for default start
2839 : //cerr<<"freqmin="<<freqmin<<" datafstart="<<datafstart<<" freqmax="<<freqmax<<" datafend="<<datafend<<endl;
2840 0 : return buildCoordinateSystemCore( msobj, spwids, fld, freqmin, freqmax, datafstart, datafend );
2841 0 : }
2842 :
2843 73 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2844 : MeasurementSet& msobj,
2845 : Vector<Int> spwids, Int fld,
2846 : Double freqmin, Double freqmax,
2847 : Double datafstart, Double datafend )
2848 : {
2849 146 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2850 :
2851 73 : CoordinateSystem csys;
2852 73 : if( csysRecord.nfields()!=0 )
2853 : {
2854 : //use cysRecord
2855 0 : Record subRec1;
2856 : // cout<<"USE THE EXISTING CSYS +++++++++++++++++"<<endl;
2857 0 : CoordinateSystem *csysptr = CoordinateSystem::restore(csysRecord,"coordsys");
2858 : //csys = *csysptr;
2859 : //CoordinateSystem csys(*csysptr);
2860 0 : csys = *csysptr;
2861 :
2862 0 : }
2863 : else {
2864 73 : MSColumns msc(msobj);
2865 73 : String telescop = msc.observation().telescopeName()(0);
2866 73 : MEpoch obsEpoch = msc.timeMeas()(0);
2867 73 : MPosition obsPosition;
2868 73 : if(!(MeasTable::Observatory(obsPosition, telescop)))
2869 : {
2870 : os << LogIO::WARN << "Did not get the position of " << telescop
2871 0 : << " from data repository" << LogIO::POST;
2872 : os << LogIO::WARN
2873 : << "Please contact CASA to add it to the repository."
2874 0 : << LogIO::POST;
2875 0 : os << LogIO::WARN << "Using first antenna position as refence " << LogIO::POST;
2876 : // unknown observatory, use first antenna
2877 0 : obsPosition=msc.antenna().positionMeas()(0);
2878 : }
2879 73 : MDirection phaseCenterToUse = phaseCenter;
2880 :
2881 73 : if( phaseCenterFieldId != -1 )
2882 : {
2883 71 : MSFieldColumns msfield(msobj.field());
2884 71 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2885 : {
2886 71 : if(trackSource){
2887 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2888 : }
2889 : else{
2890 71 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2891 : }
2892 : }
2893 : else
2894 : {
2895 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2896 : }
2897 71 : }
2898 : // Setup Phase center if it is specified only by field id.
2899 :
2900 : /////////////////// Direction Coordinates
2901 73 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2902 : // Normalize correctly
2903 73 : MVAngle ra=mvPhaseCenter.get()(0);
2904 73 : ra(0.0);
2905 :
2906 73 : MVAngle dec=mvPhaseCenter.get()(1);
2907 73 : Vector<Double> refCoord(2);
2908 73 : refCoord(0)=ra.get().getValue();
2909 73 : refCoord(1)=dec;
2910 73 : Vector<Double> refPixel(2);
2911 73 : refPixel(0) = Double(imsize[0]/2);
2912 73 : refPixel(1) = Double(imsize[1]/2);
2913 :
2914 73 : Vector<Double> deltas(2);
2915 73 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2916 73 : deltas(1)=cellsize[1].get("rad").getValue();
2917 73 : Matrix<Double> xform(2,2);
2918 73 : xform=0.0;xform.diagonal()=1.0;
2919 :
2920 73 : Vector<Double> projparams(2);
2921 73 : projparams = 0.0;
2922 73 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2923 73 : Projection projTo( projection.type(), projparams );
2924 :
2925 : DirectionCoordinate
2926 73 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2927 : // projection,
2928 : projTo,
2929 219 : refCoord(0), refCoord(1),
2930 219 : deltas(0), deltas(1),
2931 : xform,
2932 146 : refPixel(0), refPixel(1));
2933 :
2934 :
2935 : //defining observatory...needed for position on earth
2936 : // get the first ms for multiple MSes
2937 :
2938 :
2939 73 : obslocation=obsPosition;
2940 73 : ObsInfo myobsinfo;
2941 73 : myobsinfo.setTelescope(telescop);
2942 73 : myobsinfo.setPointingCenter(mvPhaseCenter);
2943 73 : myobsinfo.setObsDate(obsEpoch);
2944 73 : myobsinfo.setObserver(msc.observation().observer()(0));
2945 :
2946 : /// Attach obsInfo to the CoordinateSystem
2947 : ///csys.setObsInfo(myobsinfo);
2948 :
2949 :
2950 : /////////////////// Spectral Coordinate
2951 :
2952 : //Make sure frame conversion is switched off for REST frame data.
2953 : //Bool freqFrameValid=(freqFrame != MFrequency::REST);
2954 :
2955 : //freqFrameValid=(freqFrame != MFrequency::REST );
2956 : //UR//freqFrameValid=(freqFrame != MFrequency::REST || mode != "cubedata" );
2957 : //UR - moved freqFrameValid calc to vi/vi2 dependent wrappers.
2958 :
2959 73 : if(spwids.nelements()==0)
2960 : {
2961 0 : Int nspw=msc.spectralWindow().nrow();
2962 0 : spwids.resize(nspw);
2963 0 : indgen(spwids);
2964 : }
2965 73 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2966 73 : Vector<Double> dataChanFreq, dataChanWidth;
2967 73 : std::vector<std::vector<Int> > averageWhichChan;
2968 73 : std::vector<std::vector<Int> > averageWhichSPW;
2969 73 : std::vector<std::vector<Double> > averageChanFrac;
2970 :
2971 73 : if(spwids.nelements()==1)
2972 : {
2973 73 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2974 73 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
2975 : }
2976 : else
2977 : {
2978 0 : if(!MSTransformRegridder::combineSpwsCore(os,msobj, spwids,dataChanFreq,dataChanWidth,
2979 : averageWhichChan,averageWhichSPW,averageChanFrac))
2980 : {
2981 0 : os << LogIO::SEVERE << "Error combining SpWs" << LogIO::POST;
2982 : }
2983 : }
2984 73 : Double minDataFreq = min(dataChanFreq);
2985 73 : if(start=="" && minDataFreq < datafstart ) {
2986 : // limit data chan freq vector for default start case with channel selection
2987 : Int chanStart, chanEnd;
2988 0 : Int lochan = 0;
2989 0 : Int nDataChan = dataChanFreq.nelements();
2990 0 : Int hichan = nDataChan-1;
2991 : Double diff_fmin, diff_fmax;
2992 0 : Bool ascending = dataChanFreq[nDataChan-1] - dataChanFreq[0] > 0;
2993 0 : for(Int ichan = 0; ichan < nDataChan; ichan++)
2994 : {
2995 0 : diff_fmin = dataChanFreq[ichan] - datafstart;
2996 0 : diff_fmax = datafend - dataChanFreq[ichan];
2997 : // freqmin and freqmax should corresponds to the channel edges
2998 0 : if(ascending)
2999 : {
3000 :
3001 0 : if( diff_fmin > 0 && diff_fmin <= dataChanWidth[ichan]/2. )
3002 : {
3003 0 : lochan = ichan;
3004 : }
3005 0 : else if(diff_fmax > 0 && diff_fmax <= dataChanWidth[ichan]/2. )
3006 : {
3007 0 : hichan = ichan;
3008 : }
3009 : }
3010 : else
3011 : {
3012 0 : if( diff_fmax > 0 && diff_fmax <= dataChanWidth[ichan]/2. )
3013 : {
3014 0 : hichan = ichan;
3015 : }
3016 0 : else if( diff_fmin > 0 && diff_fmin <= dataChanWidth[ichan]/2. )
3017 : {
3018 0 : lochan = ichan;
3019 : }
3020 : }
3021 : }
3022 0 : chanStart = lochan;
3023 0 : chanEnd = hichan;
3024 0 : if (lochan > hichan)
3025 : {
3026 0 : chanStart=hichan;
3027 0 : chanEnd=lochan;
3028 : }
3029 0 : Vector<Double> tempChanFreq = dataChanFreq(Slice(chanStart,chanEnd-chanStart+1,1));
3030 0 : Vector<Double> tempChanWidth = dataChanWidth(Slice(chanStart,chanEnd-chanStart+1,1));
3031 0 : dataChanFreq.resize(tempChanFreq.nelements());
3032 0 : dataChanWidth.resize(tempChanWidth.nelements());
3033 0 : dataChanFreq = tempChanFreq;
3034 0 : dataChanWidth = tempChanWidth;
3035 0 : }
3036 73 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3037 73 : String cubemode;
3038 73 : if ( qrestfreq.getValue("Hz")==0 )
3039 : {
3040 73 : MSDopplerUtil msdoppler(msobj);
3041 73 : Vector<Double> restfreqvec;
3042 73 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3043 73 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3044 73 : if ( qrestfreq.getValue("Hz")==0 and mode!="mfs" )
3045 : {
3046 0 : cubemode = findSpecMode(mode);
3047 0 : if ( cubemode=="channel" || cubemode=="frequency" )
3048 : {
3049 : //Double provisional_restfreq = msc.spectralWindow().refFrequency()(spwids(0));
3050 : //By PLWG request, changed to center (mean) frequency of the selected spws -2015-06-22(TT)
3051 0 : Double provisional_restfreq = (datafend+datafstart)/2.0;
3052 0 : qrestfreq = Quantity(provisional_restfreq, "Hz");
3053 : os << LogIO::WARN << "No rest frequency info, using the center of the selected spw(s):"
3054 : << provisional_restfreq <<" Hz. Velocity labelling may not be correct."
3055 0 : << LogIO::POST;
3056 : }
3057 : else { // must be vel mode
3058 0 : throw(AipsError("No valid rest frequency is defined in the data, please specify the restfreq parameter") );
3059 : }
3060 : }
3061 73 : }
3062 : Double refPix;
3063 73 : Vector<Double> chanFreq;
3064 73 : Vector<Double> chanFreqStep;
3065 73 : String specmode;
3066 :
3067 73 : if(mode=="cubesource"){
3068 0 : MDoppler mdop(sysvelvalue, MDoppler::RELATIVISTIC);
3069 0 : dataChanFreq=mdop.shiftFrequency(dataChanFreq);
3070 0 : dataChanWidth=mdop.shiftFrequency(dataChanWidth);
3071 0 : if (std::isnan(dataChanFreq[0]) || std::isnan(dataChanFreq[dataChanFreq.nelements()-1])) {
3072 0 : throw(AipsError("The Doppler shift correction of the data channel frequencies resulted in 'NaN' using the radial velocity = "+
3073 0 : String::toString(sysvelvalue)+". Typically this indicates a problem in the ephemeris data being used."));
3074 : }
3075 0 : }
3076 :
3077 73 : if (!getImFreq(chanFreq, chanFreqStep, refPix, specmode, obsEpoch,
3078 : obsPosition, dataChanFreq, dataChanWidth, dataFrame, qrestfreq, freqmin, freqmax,
3079 : phaseCenterToUse))
3080 0 : throw(AipsError("Failed to determine channelization parameters"));
3081 :
3082 73 : Bool nonLinearFreq(false);
3083 73 : String veltype_p=veltype;
3084 73 : veltype_p.upcase();
3085 219 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3086 219 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3087 : {
3088 0 : nonLinearFreq= true;
3089 : }
3090 :
3091 73 : SpectralCoordinate mySpectral;
3092 : Double stepf;
3093 73 : if(!nonLinearFreq)
3094 : //TODO: velocity mode default start case (use last channels?)
3095 : {
3096 73 : Double startf=chanFreq[0];
3097 : //Double stepf=chanFreqStep[0];
3098 73 : if(chanFreq.nelements()==1)
3099 : {
3100 73 : stepf=chanFreqStep[0];
3101 : }
3102 : else
3103 : {
3104 0 : stepf=chanFreq[1]-chanFreq[0];
3105 : }
3106 73 : Double restf=qrestfreq.getValue("Hz");
3107 : //stepf=9e8;
3108 73 : if ( mode=="mfs" and restf == 0.0 ) restf = restFreq[0].getValue("Hz");
3109 : //cerr<<" startf="<<startf<<" stepf="<<stepf<<" refPix="<<refPix<<" restF="<<restf<<endl;
3110 : // once NOFRAME is implemented do this
3111 73 : if(mode=="cubedata")
3112 : {
3113 : // mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3114 0 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST?
3115 : MFrequency::REST : MFrequency::Undefined,
3116 0 : startf, stepf, refPix, restf);
3117 : }
3118 73 : else if(mode=="cubesource")
3119 : {
3120 : /*stepf=chanFreq.nelements() > 1 ?(freqmax-freqmin)/Double(chanFreq.nelements()-1) : freqmax-freqmin;
3121 : startf=freqmin+stepf/2.0;
3122 : */
3123 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
3124 0 : startf, stepf, refPix, restf);
3125 : }
3126 : else
3127 : {
3128 146 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3129 73 : startf, stepf, refPix, restf);
3130 : }
3131 : }
3132 : else
3133 : { // nonlinear freq coords - use tabular setting
3134 : // once NOFRAME is implemented do this
3135 0 : if(mode=="cubedata")
3136 : {
3137 : //mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3138 0 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST ?
3139 : MFrequency::REST : MFrequency::Undefined,
3140 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3141 : }
3142 0 : else if (mode=="cubesource")
3143 : {
3144 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
3145 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3146 : }
3147 : else
3148 : {
3149 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3150 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3151 : }
3152 : }
3153 : //cout << "Rest Freq : " << restFreq << endl;
3154 :
3155 : //for(uInt k=1 ; k < restFreq.nelements(); ++k)
3156 : //mySpectral.setRestFrequency(restFreq[k].getValue("Hz"));
3157 :
3158 73 : uInt nrestfreq = restFreq.nelements();
3159 73 : if ( nrestfreq > 1 ) {
3160 0 : Vector<Double> restfreqval( nrestfreq - 1 );
3161 0 : for ( uInt k=1 ; k < nrestfreq; ++k ) {
3162 0 : restfreqval[k-1] = restFreq[k].getValue("Hz");
3163 : }
3164 0 : mySpectral.setRestFrequencies(restfreqval, 0, true);
3165 0 : }
3166 :
3167 : // no longer needed, done inside SIImageStore
3168 : //if ( freqFrameValid ) {
3169 : // mySpectral.setReferenceConversion(MFrequency::LSRK,obsEpoch,obsPosition,phaseCenterToUse);
3170 : //}
3171 :
3172 : // cout << "RF from coordinate : " << mySpectral.restFrequency() << endl;
3173 :
3174 : ////////////////// Stokes Coordinate
3175 :
3176 73 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3177 73 : if(whichStokes.nelements()==0)
3178 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3179 73 : StokesCoordinate myStokes(whichStokes);
3180 :
3181 : ////////////////// Build Full coordinate system.
3182 :
3183 : //CoordinateSystem csys;
3184 73 : csys.addCoordinate(myRaDec);
3185 73 : csys.addCoordinate(myStokes);
3186 73 : csys.addCoordinate(mySpectral);
3187 73 : csys.setObsInfo(myobsinfo);
3188 :
3189 : //store back csys to impars record
3190 : //cerr<<"save csys to csysRecord..."<<endl;
3191 73 : if(csysRecord.isDefined("coordsys"))
3192 0 : csysRecord.removeField("coordsys");
3193 73 : csys.save(csysRecord,"coordsys");
3194 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3195 : // imshape
3196 73 : imshape.resize(4);
3197 73 : imshape[0] = imsize[0];
3198 73 : imshape[1] = imsize[0];
3199 73 : imshape[2] = whichStokes.nelements();
3200 73 : imshape[3] = chanFreq.nelements();
3201 : //toRecord();
3202 : //////////////// Set Observatory info, if MS is provided.
3203 : // (remove this section after verified...)
3204 : /***
3205 : if( ! msobj.isNull() )
3206 : {
3207 : //defining observatory...needed for position on earth
3208 : MSColumns msc(msobj);
3209 : String telescop = msc.observation().telescopeName()(0);
3210 : MEpoch obsEpoch = msc.timeMeas()(0);
3211 : MPosition obsPosition;
3212 : if(!(MeasTable::Observatory(obsPosition, telescop)))
3213 : {
3214 : os << LogIO::WARN << "Did not get the position of " << telescop
3215 : << " from data repository" << LogIO::POST;
3216 : os << LogIO::WARN
3217 : << "Please contact CASA to add it to the repository."
3218 : << LogIO::POST;
3219 : os << LogIO::WARN << "Frequency conversion will not work " << LogIO::POST;
3220 : }
3221 :
3222 : ObsInfo myobsinfo;
3223 : myobsinfo.setTelescope(telescop);
3224 : myobsinfo.setPointingCenter(mvPhaseCenter);
3225 : myobsinfo.setObsDate(obsEpoch);
3226 : myobsinfo.setObserver(msc.observation().observer()(0));
3227 :
3228 : /// Attach obsInfo to the CoordinateSystem
3229 : csys.setObsInfo(myobsinfo);
3230 :
3231 : }// if MS is provided.
3232 : ***/
3233 73 : } // end of else when coordsys record is not defined...
3234 :
3235 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3236 :
3237 146 : return csys;
3238 73 : }
3239 :
3240 :
3241 : /*
3242 : #ifdef USEVIVB2
3243 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystem(vi::VisibilityIterator2* vi2)
3244 : #else
3245 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystem(ROVisibilityIterator* rvi )
3246 : #endif
3247 : {
3248 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
3249 :
3250 :
3251 : // get the first ms for multiple MSes
3252 : #ifdef USEVIVB2
3253 : MeasurementSet msobj=vi2->getMeasurementSet();
3254 : #else
3255 : MeasurementSet msobj=rvi->getMeasurementSet();
3256 : #endif
3257 :
3258 : MDirection phaseCenterToUse = phaseCenter;
3259 : if( phaseCenterFieldId != -1 )
3260 : {
3261 : MSFieldColumns msfield(msobj.field());
3262 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
3263 : }
3264 : // Setup Phase center if it is specified only by field id.
3265 :
3266 : /////////////////// Direction Coordinates
3267 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
3268 : // Normalize correctly
3269 : MVAngle ra=mvPhaseCenter.get()(0);
3270 : ra(0.0);
3271 :
3272 : MVAngle dec=mvPhaseCenter.get()(1);
3273 : Vector<Double> refCoord(2);
3274 : refCoord(0)=ra.get().getValue();
3275 : refCoord(1)=dec;
3276 : Vector<Double> refPixel(2);
3277 : refPixel(0) = Double(imsize[0]/2);
3278 : refPixel(1) = Double(imsize[1]/2);
3279 :
3280 : Vector<Double> deltas(2);
3281 : deltas(0)=-1* cellsize[0].get("rad").getValue();
3282 : deltas(1)=cellsize[1].get("rad").getValue();
3283 : Matrix<Double> xform(2,2);
3284 : xform=0.0;xform.diagonal()=1.0;
3285 :
3286 : Vector<Double> projparams(2);
3287 : projparams = 0.0;
3288 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
3289 : Projection projTo( projection.type(), projparams );
3290 :
3291 : DirectionCoordinate
3292 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
3293 : // projection,
3294 : projTo,
3295 : refCoord(0), refCoord(1),
3296 : deltas(0), deltas(1),
3297 : xform,
3298 : refPixel(0), refPixel(1));
3299 :
3300 :
3301 : //defining observatory...needed for position on earth
3302 : // get the first ms for multiple MSes
3303 : MSColumns msc(msobj);
3304 : String telescop = msc.observation().telescopeName()(0);
3305 : MEpoch obsEpoch = msc.timeMeas()(0);
3306 : MPosition obsPosition;
3307 : if(!(MeasTable::Observatory(obsPosition, telescop)))
3308 : {
3309 : os << LogIO::WARN << "Did not get the position of " << telescop
3310 : << " from data repository" << LogIO::POST;
3311 : os << LogIO::WARN
3312 : << "Please contact CASA to add it to the repository."
3313 : << LogIO::POST;
3314 : os << LogIO::WARN << "Frequency conversion will not work " << LogIO::POST;
3315 : }
3316 :
3317 : ObsInfo myobsinfo;
3318 : myobsinfo.setTelescope(telescop);
3319 : myobsinfo.setPointingCenter(mvPhaseCenter);
3320 : myobsinfo.setObsDate(obsEpoch);
3321 : myobsinfo.setObserver(msc.observation().observer()(0));
3322 :
3323 : /// Attach obsInfo to the CoordinateSystem
3324 : ///csys.setObsInfo(myobsinfo);
3325 :
3326 :
3327 : /////////////////// Spectral Coordinate
3328 :
3329 : //Make sure frame conversion is switched off for REST frame data.
3330 : //Bool freqFrameValid=(freqFrame != MFrequency::REST);
3331 :
3332 : //freqFrameValid=(freqFrame != MFrequency::REST );
3333 : freqFrameValid=(freqFrame != MFrequency::REST || mode != "cubedata" );
3334 :
3335 : // *** get selected spw ids ***
3336 : Vector<Int> spwids;
3337 : #ifdef USEVIVB2
3338 : vi2->spectralWindows( spwids );
3339 : #else
3340 : Vector<Int> nvischan;
3341 : rvi->allSelectedSpectralWindows(spwids,nvischan);
3342 : #endif
3343 : if(spwids.nelements()==0)
3344 : {
3345 : Int nspw=msc.spectralWindow().nrow();
3346 : spwids.resize(nspw);
3347 : indgen(spwids);
3348 : }
3349 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
3350 : Vector<Double> dataChanFreq, dataChanWidth;
3351 : if(spwids.nelements()==1)
3352 : {
3353 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
3354 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
3355 : }
3356 : else
3357 : {
3358 : SubMS thems(msobj);
3359 : if(!thems.combineSpws(spwids,true,dataChanFreq,dataChanWidth))
3360 : //if(!MSTransformRegridder::combineSpws(os,msobj.tableName(),spwids,dataChanFreq,dataChanWidth))
3361 : {
3362 : os << LogIO::SEVERE << "Error combining SpWs" << LogIO::POST;
3363 : }
3364 : }
3365 :
3366 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3367 : if( qrestfreq.getValue("Hz")==0 )
3368 : {
3369 : #ifdef USEVIVB2
3370 : ///// TTCheck
3371 : Vector<Int> flds;
3372 : vi2->fieldIds( flds );
3373 : AlwaysAssert( flds.nelements()>0 , AipsError );
3374 : Int fld = flds[0];
3375 : #else
3376 : Int fld = rvi->fieldId();
3377 : #endif
3378 : MSDopplerUtil msdoppler(msobj);
3379 : Vector<Double> restfreqvec;
3380 : msdoppler.dopplerInfo(restfreqvec, spwids[0], fld);
3381 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec[0],"Hz"): Quantity(0.0, "Hz");
3382 : }
3383 : Double refPix;
3384 : Vector<Double> chanFreq;
3385 : Vector<Double> chanFreqStep;
3386 : String specmode;
3387 :
3388 : //for mfs
3389 : Double freqmin=0, freqmax=0;
3390 : #ifdef USEVIVB2
3391 : vi2->getFreqInSpwRange(freqmin,freqmax,freqFrameValid? freqFrame:MFrequency::REST );
3392 : #else
3393 : rvi->getFreqInSpwRange(freqmin,freqmax,freqFrameValid? freqFrame:MFrequency::REST );
3394 : #endif
3395 :
3396 : if (!getImFreq(chanFreq, chanFreqStep, refPix, specmode, obsEpoch,
3397 : obsPosition, dataChanFreq, dataChanWidth, dataFrame, qrestfreq, freqmin, freqmax,
3398 : phaseCenterToUse))
3399 : throw(AipsError("Failed to determine channelization parameters"));
3400 :
3401 : Bool nonLinearFreq(false);
3402 : String veltype_p=veltype;
3403 : veltype_p.upcase();
3404 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3405 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3406 : {
3407 : nonLinearFreq= true;
3408 : }
3409 :
3410 : SpectralCoordinate mySpectral;
3411 : Double stepf;
3412 : if(!nonLinearFreq)
3413 : //TODO: velocity mode default start case (use last channels?)
3414 : {
3415 : Double startf=chanFreq[0];
3416 : //Double stepf=chanFreqStep[0];
3417 : if(chanFreq.nelements()==1)
3418 : {
3419 : stepf=chanFreqStep[0];
3420 : }
3421 : else
3422 : {
3423 : stepf=chanFreq[1]-chanFreq[0];
3424 : }
3425 : Double restf=qrestfreq.getValue("Hz");
3426 : //cerr<<" startf="<<startf<<" stepf="<<stepf<<" refPix="<<refPix<<" restF="<<restf<<endl;
3427 : // once NOFRAME is implemented do this
3428 : if(mode=="cubedata")
3429 : {
3430 : // mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3431 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST?
3432 : MFrequency::REST : MFrequency::Undefined,
3433 : startf, stepf, refPix, restf);
3434 : }
3435 : else
3436 : {
3437 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3438 : startf, stepf, refPix, restf);
3439 : }
3440 : }
3441 : else
3442 : { // nonlinear freq coords - use tabular setting
3443 : // once NOFRAME is implemented do this
3444 : if(mode=="cubedata")
3445 : {
3446 : //mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3447 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST ?
3448 : MFrequency::REST : MFrequency::Undefined,
3449 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3450 : }
3451 : else
3452 : {
3453 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3454 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3455 : }
3456 : }
3457 : //cout << "Rest Freq : " << restFreq << endl;
3458 :
3459 : for(uInt k=1 ; k < restFreq.nelements(); ++k)
3460 : mySpectral.setRestFrequency(restFreq[k].getValue("Hz"));
3461 :
3462 : if ( freqFrameValid ) {
3463 : mySpectral.setReferenceConversion(MFrequency::LSRK,obsEpoch,obsPosition,phaseCenterToUse);
3464 : }
3465 :
3466 : // cout << "RF from coordinate : " << mySpectral.restFrequency() << endl;
3467 :
3468 : ////////////////// Stokes Coordinate
3469 :
3470 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3471 : if(whichStokes.nelements()==0)
3472 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3473 : StokesCoordinate myStokes(whichStokes);
3474 :
3475 : ////////////////// Build Full coordinate system.
3476 :
3477 : CoordinateSystem csys;
3478 : csys.addCoordinate(myRaDec);
3479 : csys.addCoordinate(myStokes);
3480 : csys.addCoordinate(mySpectral);
3481 : csys.setObsInfo(myobsinfo);
3482 :
3483 : //////////////// Set Observatory info, if MS is provided.
3484 : // (remove this section after verified...)
3485 : return csys;
3486 : }
3487 : */
3488 :
3489 0 : MDirection SynthesisParamsImage::getMovingSourceDir(const MeasurementSet& ms, const MEpoch& refEp, const MPosition& obsposition, const MDirection::Types outframe){
3490 0 : MDirection outdir;
3491 0 : String ephemtab(movingSource);
3492 0 : if(movingSource=="TRACKFIELD"){
3493 0 : Int fieldID=MSColumns(ms).fieldId()(0);
3494 0 : ephemtab=Path(MSColumns(ms).field().ephemPath(fieldID)).absoluteName();
3495 : }
3496 0 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
3497 0 : if( (! Table::isReadable(ephemtab)) && ( (planetType <= MDirection::N_Types) || (planetType >= MDirection::COMET)))
3498 0 : throw(AipsError("Does not have a valid ephemeris table or major solar system object defined"));
3499 0 : MeasFrame mframe(refEp, obsposition);
3500 0 : MDirection::Ref outref1(MDirection::AZEL, mframe);
3501 0 : MDirection::Ref outref(outframe, mframe);
3502 0 : MDirection tmpazel;
3503 : // (TT) Switched the order of evaluation of if statement (if ephem table is readable
3504 : // one should use that. MDirection::getType will match MDirection::Types if a string conains and starts with
3505 : // one of the enum names in MDirection::Types. So the table name can be mistaken as a major planets in MDirection::Types
3506 : // if it is evaluated first.
3507 0 : if (Table::isReadable(ephemtab)){
3508 0 : MeasComet mcomet(Path(ephemtab).absoluteName());
3509 0 : mframe.set(mcomet);
3510 0 : tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
3511 0 : }
3512 0 : else if (planetType >=MDirection::MERCURY && planetType <MDirection::COMET){
3513 0 : tmpazel=MDirection::Convert(trackDir, outref1)();
3514 : }
3515 0 : outdir=MDirection::Convert(tmpazel, outref)();
3516 :
3517 0 : return outdir;
3518 0 : }
3519 :
3520 73 : Bool SynthesisParamsImage::getImFreq(Vector<Double>& chanFreq, Vector<Double>& chanFreqStep,
3521 : Double& refPix, String& specmode,
3522 : const MEpoch& obsEpoch, const MPosition& obsPosition,
3523 : const Vector<Double>& dataChanFreq,
3524 : const Vector<Double>& dataChanWidth,
3525 : const MFrequency::Types& dataFrame,
3526 : const Quantity& qrestfreq, const Double& freqmin, const Double& freqmax,
3527 : const MDirection& phaseCenter)
3528 : {
3529 :
3530 73 : String inStart, inStep;
3531 73 : specmode = findSpecMode(mode);
3532 73 : String freqframe;
3533 73 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3534 146 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3535 :
3536 73 : refPix=0.0;
3537 73 : Bool descendingfreq(false);
3538 73 : Bool descendingoutfreq(false);
3539 :
3540 73 : if( mode.contains("cube") )
3541 : {
3542 0 : String restfreq=QuantityToString(qrestfreq);
3543 : // use frame from input start or width in MFreaquency or MRadialVelocity
3544 0 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3545 : // emit warning here if qmframe is used
3546 : //
3547 0 : inStart = start;
3548 0 : inStep = step;
3549 0 : if( specmode=="channel" )
3550 : {
3551 0 : inStart = String::toString(chanStart);
3552 0 : inStep = String::toString(chanStep);
3553 : // negative step -> descending channel indices
3554 0 : if (inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3555 : // input frame is the data frame
3556 : //freqframe = MFrequency::showType(dataFrame);
3557 : }
3558 0 : else if( specmode=="frequency" )
3559 : {
3560 : //if ( freqStart.getValue("Hz") == 0 && freqStart.getUnit() != "" ) { // default start
3561 : //start = String::toString( freqmin ) + freqStart.getUnit();
3562 : //}
3563 : //else {
3564 : //start = String::toString( freqStart.getValue(freqStart.getUnit()) )+freqStart.getUnit();
3565 : //}
3566 : //step = String::toString( freqStep.getValue(freqStep.getUnit()) )+freqStep.getUnit();
3567 : // negative freq width -> descending freq ordering
3568 0 : if(inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3569 : }
3570 0 : else if( specmode=="velocity" )
3571 : {
3572 : // if velStart is empty set start to vel of freqmin or freqmax?
3573 : //if ( velStart.getValue(velStart.getUnit()) == 0 && !(velStart.getUnit().contains("m/s")) ) {
3574 : // start = "";
3575 : //}
3576 : //else {
3577 : // start = String::toString( velStart.getValue(velStart.getUnit()) )+velStart.getUnit();
3578 : //}
3579 : //step = String::toString( velStep.getValue(velStep.getUnit()) )+velStep.getUnit();
3580 : // positive velocity width -> descending freq ordering
3581 0 : if (!inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3582 : }
3583 :
3584 0 : if (inStep=='0') inStep="";
3585 :
3586 0 : MRadialVelocity mSysVel;
3587 0 : Quantity qVel;
3588 : MRadialVelocity::Types mRef;
3589 0 : if(mode!="cubesource")
3590 : {
3591 :
3592 :
3593 0 : if(freqframe=="SOURCE")
3594 : {
3595 : os << LogIO::SEVERE << "freqframe=\"SOURCE\" is only allowed for mode=\"cubesrc\""
3596 0 : << LogIO::EXCEPTION;
3597 0 : return false;
3598 : }
3599 : }
3600 : else // only for cubesrc mode: TODO- check for the ephemeris info.
3601 : {
3602 0 : freqframe=MFrequency::showType(dataFrame);
3603 0 : if(sysvel!="") {
3604 0 : stringToQuantity(sysvel,qVel);
3605 0 : MRadialVelocity::getType(mRef,sysvelframe);
3606 0 : mSysVel=MRadialVelocity(qVel,mRef);
3607 : }
3608 : else // and if no ephemeris info, issue a warning...
3609 0 : { mSysVel=MRadialVelocity();}
3610 : }
3611 : // cubedata mode: input start, step are those of the input data frame
3612 0 : if ( mode=="cubedata" )
3613 : {
3614 0 : freqframe=MFrequency::showType(dataFrame);
3615 0 : freqFrameValid=false; // no conversion for vb.lsrfrequency()
3616 : }
3617 : //if ( mode=="cubedata" ) freqframe=MFrequency::REST;
3618 :
3619 : // *** NOTE ***
3620 : // calcChanFreqs alway returns chanFreq in
3621 : // ascending freq order.
3622 : // for step < 0 calcChanFreqs returns chanFreq that
3623 : // contains start freq. in its last element of the vector.
3624 : //
3625 0 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3626 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3627 0 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3628 0 : << LogIO::POST;
3629 0 : ostringstream ostr;
3630 0 : ostr << " phaseCenter='" << phaseCenter;
3631 0 : os << String(ostr)<<"' ";
3632 :
3633 : Double dummy; // dummy variable - weightScale is not used here
3634 0 : Bool rst=MSTransformRegridder::calcChanFreqs(os,
3635 : chanFreq,
3636 : chanFreqStep,
3637 : dummy,
3638 : dataChanFreq,
3639 : dataChanWidth,
3640 : phaseCenter,
3641 : dataFrame,
3642 : obsEpoch,
3643 : obsPosition,
3644 : specmode,
3645 : nchan,
3646 : inStart,
3647 : inStep,
3648 : restfreq,
3649 : freqframe,
3650 0 : veltype,
3651 : verbose,
3652 : mSysVel
3653 : );
3654 :
3655 0 : if( nchan==-1 )
3656 : {
3657 0 : nchan = chanFreq.nelements();
3658 0 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3659 : }
3660 :
3661 0 : if (!rst) {
3662 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
3663 0 : << LogIO::EXCEPTION;
3664 0 : return false;
3665 : }
3666 : os << LogIO::DEBUG1
3667 0 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3668 0 : << LogIO::POST;
3669 :
3670 0 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3671 0 : descendingoutfreq = true;
3672 : }
3673 :
3674 : //if (descendingfreq && !descendingoutfreq) {
3675 0 : if ((specmode=="channel" && descendingfreq==1)
3676 0 : || (specmode!="channel" && (descendingfreq != descendingoutfreq))) {
3677 : // reverse the freq vector if necessary so the first element can be
3678 : // used to set spectralCoordinates in all the cases.
3679 : //
3680 : // also do for chanFreqStep..
3681 0 : std::vector<Double> stlchanfreq;
3682 0 : chanFreq.tovector(stlchanfreq);
3683 0 : std::reverse(stlchanfreq.begin(),stlchanfreq.end());
3684 0 : chanFreq=Vector<Double>(stlchanfreq);
3685 0 : chanFreqStep=-chanFreqStep;
3686 0 : }
3687 0 : }
3688 73 : else if ( mode=="mfs" ) {
3689 73 : chanFreq.resize(1);
3690 73 : chanFreqStep.resize(1);
3691 : //chanFreqStep[0] = freqmax - freqmin;
3692 73 : Double freqmean = (freqmin + freqmax)/2;
3693 73 : if (refFreq.getValue("Hz")==0) {
3694 73 : chanFreq[0] = freqmean;
3695 73 : refPix = 0.0;
3696 73 : chanFreqStep[0] = freqmax - freqmin;
3697 : }
3698 : else {
3699 0 : chanFreq[0] = refFreq.getValue("Hz");
3700 : // Set the new reffreq to be the refPix (CAS-9518)
3701 0 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3702 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3703 0 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3704 : }
3705 :
3706 73 : if( nchan==-1 ) nchan=1;
3707 73 : if( qrestfreq.getValue("Hz")==0.0 ) {
3708 0 : restFreq.resize(1);
3709 0 : restFreq[0] = Quantity(freqmean,"Hz");
3710 : }
3711 : }
3712 : else {
3713 : // unrecognized mode, error
3714 0 : os << LogIO::SEVERE << "mode="<<mode<<" is unrecognized."
3715 0 : << LogIO::EXCEPTION;
3716 0 : return false;
3717 : }
3718 73 : return true;
3719 :
3720 73 : }//getImFreq
3721 : /////////////////////////
3722 73 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3723 73 : Double inStartFreq=-1.0;
3724 73 : String checkspecmode("");
3725 73 : if(mode.contains("cube")) {
3726 0 : checkspecmode = findSpecMode(mode);
3727 : }
3728 73 : if(checkspecmode!="") {
3729 0 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3730 0 : if(checkspecmode=="channel") {
3731 0 : inStartFreq=-1.0;
3732 : }
3733 : else {
3734 0 : if(checkspecmode=="frequency") {
3735 0 : inStartFreq = freqStart.get("Hz").getValue();
3736 : }
3737 0 : else if(checkspecmode=="velocity") {
3738 : MDoppler::Types DopType;
3739 0 : MDoppler::getType(DopType, veltype);
3740 0 : MDoppler mdop(velStart,DopType);
3741 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3742 0 : inStartFreq = MFrequency::fromDoppler(mdop, qrestfreq.getValue(Unit("Hz")), mfreqframe).getValue();
3743 0 : }
3744 : }
3745 : }
3746 :
3747 73 : return inStartFreq;
3748 :
3749 73 : }
3750 :
3751 73 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3752 : {
3753 73 : String specmode;
3754 73 : specmode="channel";
3755 73 : if ( mode.contains("cube") ) {
3756 : // if velstart or velstep is defined -> specmode='vel'
3757 : // else if freqstart or freqstep is defined -> specmode='freq'
3758 : // velocity: assume unset if velStart => 0.0 with no unit
3759 : // assume unset if velStep => 0.0 with/without unit
3760 0 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3761 0 : !( velStep.getValue()==0.0)) {
3762 0 : specmode="velocity";
3763 : }
3764 0 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3765 0 : !(freqStep.getValue()==0.0)) {
3766 0 : specmode="frequency";
3767 : }
3768 : }
3769 73 : return specmode;
3770 0 : }
3771 :
3772 :
3773 219 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3774 : {
3775 219 : Vector<Int> whichStokes(0);
3776 219 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3777 0 : stokes=="RR" ||stokes=="LL" ||
3778 219 : stokes=="XX" || stokes=="YY" ) {
3779 219 : whichStokes.resize(1);
3780 219 : whichStokes(0)=Stokes::type(stokes);
3781 : }
3782 0 : else if(stokes=="IV" || stokes=="IQ" ||
3783 0 : stokes=="RRLL" || stokes=="XXYY" ||
3784 0 : stokes=="QU" || stokes=="UV"){
3785 0 : whichStokes.resize(2);
3786 :
3787 0 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3788 0 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3789 0 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3790 0 : else if(stokes=="XXYY"){whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
3791 0 : else if(stokes=="QU"){whichStokes[0]=Stokes::Q; whichStokes[1]=Stokes::U; }
3792 0 : else if(stokes=="UV"){ whichStokes[0]=Stokes::U; whichStokes[1]=Stokes::V; }
3793 :
3794 : }
3795 :
3796 0 : else if(stokes=="IQU" || stokes=="IUV") {
3797 0 : whichStokes.resize(3);
3798 0 : if(stokes=="IUV")
3799 0 : {whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::U; whichStokes[2]=Stokes::V;}
3800 : else
3801 0 : {whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q; whichStokes[2]=Stokes::U;}
3802 : }
3803 0 : else if(stokes=="IQUV"){
3804 0 : whichStokes.resize(4);
3805 0 : whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;
3806 0 : whichStokes(2)=Stokes::U; whichStokes(3)=Stokes::V;
3807 : }
3808 :
3809 219 : return whichStokes;
3810 0 : }// decidenpolplanes
3811 :
3812 146 : IPosition SynthesisParamsImage::shp() const
3813 : {
3814 146 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3815 :
3816 146 : if( imsize[0]<=0 || imsize[1]<=0 || nStokes<=0 || nchan<=0 )
3817 : {
3818 0 : throw(AipsError("Internal Error : Image shape is invalid : [" + String::toString(imsize[0]) + "," + String::toString(imsize[1]) + "," + String::toString(nStokes) + "," + String::toString(nchan) + "]" ));
3819 : }
3820 :
3821 292 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3822 : }
3823 :
3824 73 : Record SynthesisParamsImage::getcsys() const
3825 : {
3826 73 : return csysRecord;
3827 : }
3828 :
3829 0 : Record SynthesisParamsImage::updateParams(const Record& impar)
3830 : {
3831 0 : Record newimpar( impar );
3832 0 : if ( impar.isDefined("csys") )
3833 : {
3834 0 : Vector<Int> newimsize(2);
3835 0 : newimsize[0] = imshape[0];
3836 0 : newimsize[1] = imshape[1];
3837 0 : newimpar.define("imsize", newimsize);
3838 0 : if ( newimpar.isDefined("direction0") )
3839 : {
3840 0 : Record dirRec = newimpar.subRecord("direction0");
3841 0 : Vector<Double> cdelta = dirRec.asArrayDouble("cdelt");
3842 0 : Vector<String> cells(2);
3843 0 : cells[0] = String::toString(-1*cdelta[0]) + "rad";
3844 0 : cells[1] = String::toString(-1*cdelta[1]) + "rad";
3845 0 : newimpar.define("cell", cells );
3846 0 : }
3847 0 : if ( newimpar.isDefined("stokes1") )
3848 : {
3849 0 : Record stokesRec = newimpar.subRecord("stokes1");
3850 0 : Vector<String> stokesvecs = stokesRec.asArrayString("stokes");
3851 0 : String stokesStr;
3852 0 : for (uInt j = 0; j < stokesvecs.nelements(); j++)
3853 : {
3854 0 : stokesStr+=stokesvecs[j];
3855 : }
3856 0 : }
3857 0 : if ( newimpar.isDefined("nchan") )
3858 : {
3859 0 : newimpar.define("nchan",imshape[2]);
3860 : }
3861 0 : if ( newimpar.isDefined("spectral2") )
3862 : {
3863 0 : Record specRec = newimpar.subRecord("spectral2");
3864 0 : if ( specRec.isDefined("restfreqs") )
3865 : {
3866 0 : Vector<Double> restfs = specRec.asArrayDouble("restfreqs");
3867 0 : Vector<String> restfstrs(restfs.nelements());
3868 0 : for(uInt restf=0; restf<restfs.nelements(); restf++){restfstrs[restf] = String::toString(restfs[restf]) + "Hz";}
3869 0 : newimpar.define("restfreq",restfstrs);
3870 0 : }
3871 : //reffreq?
3872 : //outframe
3873 : //sysvel
3874 : //sysvelframe
3875 0 : }
3876 0 : }
3877 0 : return newimpar;
3878 0 : }
3879 :
3880 : /////////////////////// Grid/FTMachine Parameters
3881 :
3882 219 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3883 : {
3884 219 : setDefaults();
3885 219 : }
3886 :
3887 219 : SynthesisParamsGrid::~SynthesisParamsGrid()
3888 : {
3889 219 : }
3890 :
3891 :
3892 73 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3893 : {
3894 73 : setDefaults();
3895 :
3896 73 : String err("");
3897 :
3898 : try {
3899 73 : err += readVal( inrec, String("imagename"), imageName );
3900 :
3901 : // FTMachine parameters
3902 73 : err += readVal( inrec, String("gridder"), gridder );
3903 73 : err += readVal( inrec, String("padding"), padding );
3904 73 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3905 73 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3906 73 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3907 73 : err += readVal( inrec, String("convfunc"), convFunc );
3908 :
3909 73 : err += readVal( inrec, String("vptable"), vpTable );
3910 :
3911 :
3912 :
3913 : // convert 'gridder' to 'ftmachine' and 'mtype'
3914 73 : ftmachine = "gridft";
3915 73 : mType = "default";
3916 73 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
3917 69 : ftmachine = "gridft";
3918 : }
3919 8 : else if ( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) &&
3920 4 : (wprojplanes>1 || wprojplanes==-1) ) {
3921 4 : ftmachine = "wprojectft";
3922 : }
3923 : //facetting alone use gridft
3924 0 : else if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes==1))
3925 0 : {ftmachine=="gridft";}
3926 :
3927 0 : else if (gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" ) {
3928 0 : ftmachine = "mosaicft";
3929 : }
3930 :
3931 0 : else if (gridder=="imagemosaic") {
3932 0 : mType = "imagemosaic";
3933 0 : if (wprojplanes>1 || wprojplanes==-1) {
3934 0 : ftmachine = "wprojectft";
3935 : }
3936 : }
3937 :
3938 0 : else if (gridder=="awproject" || gridder=="awprojectft" || gridder=="awp") {
3939 0 : ftmachine = "awprojectft";
3940 : }
3941 :
3942 0 : else if (gridder=="singledish") {
3943 :
3944 0 : ftmachine = "sd";
3945 : }
3946 : else{
3947 0 : ftmachine=gridder;
3948 0 : ftmachine.downcase();
3949 :
3950 : }
3951 :
3952 73 : String deconvolver;
3953 73 : err += readVal( inrec, String("deconvolver"), deconvolver );
3954 73 : if (deconvolver=="mtmfs") {
3955 0 : mType = "multiterm"; // Takes precedence over imagemosaic
3956 : }
3957 :
3958 : // facets
3959 73 : err += readVal( inrec, String("facets"), facets );
3960 : // chanchunks
3961 73 : err += readVal( inrec, String("chanchunks"), chanchunks );
3962 :
3963 : // Spectral interpolation
3964 73 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
3965 : // Track moving source ?
3966 73 : err += readVal( inrec, String("distance"), distance );
3967 73 : err += readVal( inrec, String("tracksource"), trackSource );
3968 73 : err += readVal( inrec, String("trackdir"), trackDir );
3969 :
3970 : // The extra params for WB-AWP
3971 73 : err += readVal( inrec, String("aterm"), aTermOn );
3972 73 : err += readVal( inrec, String("psterm"), psTermOn );
3973 73 : err += readVal( inrec, String("mterm"), mTermOn );
3974 73 : err += readVal( inrec, String("wbawp"), wbAWP );
3975 73 : err += readVal( inrec, String("cfcache"), cfCache );
3976 73 : err += readVal( inrec, String("usepointing"), usePointing );
3977 73 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3978 73 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3979 73 : err += readVal( inrec, String("conjbeams"), conjBeams );
3980 73 : err += readVal( inrec, String("computepastep"), computePAStep );
3981 73 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3982 :
3983 : // The extra params for single-dish
3984 73 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3985 73 : err += readVal( inrec, String("convertfirst"), convertFirst );
3986 73 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3987 73 : err += readVal( inrec, String("convsupport"), convSupport );
3988 73 : err += readVal( inrec, String("truncate"), truncateSize );
3989 73 : err += readVal( inrec, String("gwidth"), gwidth );
3990 73 : err += readVal( inrec, String("jwidth"), jwidth );
3991 73 : err += readVal( inrec, String("minweight"), minWeight );
3992 73 : err += readVal( inrec, String("clipminmax"), clipMinMax );
3993 :
3994 : // Single or MultiTerm mapper : read in 'deconvolver' and set mType here.
3995 : // err += readVal( inrec, String("mtype"), mType );
3996 :
3997 73 : if (ftmachine=="awprojectft" && cfCache=="") {
3998 0 : cfCache = imageName + ".cf";
3999 : }
4000 :
4001 73 : if ( ftmachine=="awprojectft" &&
4002 73 : usePointing==True &&
4003 0 : pointingOffsetSigDev.nelements() != 2 ) {
4004 : // Set the default to a large value so that it behaves like CASA 5.6's usepointing=True.
4005 0 : pointingOffsetSigDev.resize(2);
4006 0 : pointingOffsetSigDev[0] = 600.0;
4007 0 : pointingOffsetSigDev[1] = 600.0;
4008 : }
4009 :
4010 73 : err += verify();
4011 :
4012 73 : } catch(AipsError &x) {
4013 0 : err = err + x.getMesg() + "\n";
4014 0 : }
4015 :
4016 73 : if (err.length()>0) {
4017 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4018 : }
4019 :
4020 73 : }
4021 :
4022 73 : String SynthesisParamsGrid::verify() const
4023 : {
4024 73 : String err;
4025 :
4026 : // Check for valid FTMachine type.
4027 : // Valid other params per FTM type, etc... ( check about nterms>1 )
4028 :
4029 :
4030 :
4031 73 : if ( imageName == "" ) {
4032 0 : err += "Please supply an image name\n";
4033 : }
4034 :
4035 4 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4036 73 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4037 77 : (ftmachine != "mawprojectft") &&
4038 0 : (ftmachine != "sd"))
4039 : {
4040 : err += "Invalid ftmachine name. Must be one of"
4041 : " 'gridft', 'wprojectft',"
4042 : " 'mosaicft', 'awprojectft',"
4043 : " 'mawpojectft', 'protoft',"
4044 0 : " 'sd'\n";
4045 : }
4046 :
4047 :
4048 146 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4049 73 : ( ftmachine == "awprojectft" and mType == "imagemosaic" ) ) {
4050 0 : err += "Cannot use " + ftmachine + " with " + mType +
4051 : " because it is a redundant choice for mosaicing."
4052 : " In the future, we may support the combination"
4053 : " to signal the use of single-pointing sized image grids"
4054 : " during gridding and iFT,"
4055 : " and only accumulating it on the large mosaic image."
4056 : " For now, please set"
4057 : " either mappertype='default' to get mosaic gridding"
4058 0 : " or ftmachine='ft' or 'wprojectft' to get image domain mosaics.\n";
4059 : }
4060 :
4061 73 : if ( facets < 1 ) {
4062 0 : err += "Must have at least 1 facet\n";
4063 : }
4064 :
4065 : //if( chanchunks < 1 )
4066 : // {err += "Must have at least 1 chanchunk\n"; }
4067 73 : if ( facets > 1 and chanchunks > 1 ) {
4068 : err += "The combination of facetted imaging"
4069 : " with channel chunking is not yet supported."
4070 0 : " Please choose only one or the other for now.\n";
4071 : }
4072 :
4073 73 : if ( ftmachine == "wproject" and ( wprojplanes == 0 or wprojplanes == 1 ) ) {
4074 : err += "The wproject gridder must be accompanied with"
4075 0 : " wprojplanes>1 or wprojplanes=-1\n";
4076 : }
4077 :
4078 73 : if ( ftmachine == "awprojectft" and facets > 1 ) {
4079 : err += "The awprojectft gridder supports A- and W-Projection."
4080 : " Instead of using facets>1 to deal with the W-term,"
4081 : " please set the number of wprojplanes to a value > 1"
4082 0 : " to trigger the combined AW-Projection algorithm. \n";
4083 : // Also, the way the AWP cfcache is managed,
4084 : // even if all facets share a common one so that they reuse convolution functions,
4085 : // the first facet's gridder writes out the avgPB
4086 : // and all others see that it's there and don't compute their own.
4087 : // As a result, the code will run,
4088 : // but the first facet's weight image will be duplicated for all facets.
4089 : // If needed, this must be fixed in the way the AWP gridder manages its cfcache.
4090 : // But, since the AWP gridder supports joint A and W projection,
4091 : // facet support may never be needed in the first place...
4092 : }
4093 :
4094 73 : if ( ftmachine == "awprojectft" and wprojplanes == -1 ) {
4095 : err += "The awprojectft gridder does not support wprojplanes=-1"
4096 0 : " for automatic calculation. Please pick a value >1\n";
4097 : }
4098 :
4099 73 : if ( ftmachine == "mosaicft" and facets > 1 ) {
4100 : err += "The combination of mosaicft gridding"
4101 : " with multiple facets is not supported."
4102 : " Please use the awprojectft gridder instead,"
4103 0 : " and set wprojplanes to a value > 1 to trigger AW-Projection.\n";
4104 : }
4105 :
4106 73 : if ( ftmachine == "awprojectft" and usePointing == True and
4107 0 : pointingOffsetSigDev.nelements() != 2 ) {
4108 : err += "The pointingoffsetsigdev parameter must be"
4109 : " a two-element vector of doubles in order to be used with usepointing=True"
4110 0 : " and the AWProject gridder. Setting it to the default of \n";
4111 : }
4112 :
4113 : // Single-dish parameters check
4114 73 : if ( ftmachine == "sd" ) {
4115 0 : if ( convertFirst != "always" and
4116 0 : convertFirst != "never" and
4117 0 : convertFirst != "auto" ) {
4118 0 : err += "convertfirst parameter: illegal value: '" + convertFirst + "'."
4119 0 : " Allowed values: 'always', 'never', 'auto'.\n";
4120 : }
4121 : }
4122 :
4123 73 : return err;
4124 0 : }
4125 :
4126 292 : void SynthesisParamsGrid::setDefaults()
4127 : {
4128 292 : imageName="";
4129 : // FTMachine parameters
4130 : //ftmachine="GridFT";
4131 292 : ftmachine="gridft";
4132 292 : gridder=ftmachine;
4133 292 : padding=1.2;
4134 292 : useAutoCorr=false;
4135 292 : useDoublePrec=true;
4136 292 : wprojplanes=1;
4137 292 : convFunc="SF";
4138 292 : vpTable="";
4139 :
4140 : // facets
4141 292 : facets=1;
4142 :
4143 : // chanchunks
4144 292 : chanchunks=1;
4145 :
4146 : // Spectral Axis interpolation
4147 292 : interpolation=String("nearest");
4148 :
4149 : //mosaic use pointing
4150 292 : usePointing=false;
4151 : // Moving phase center ?
4152 292 : distance=Quantity(0,"m");
4153 292 : trackSource=false;
4154 292 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4155 :
4156 : // The extra params for WB-AWP
4157 292 : aTermOn = true;
4158 292 : psTermOn = true;
4159 292 : mTermOn = false;
4160 292 : wbAWP = true;
4161 292 : cfCache = "";
4162 292 : usePointing = false;
4163 292 : pointingOffsetSigDev.resize(0);
4164 : // pointingOffsetSigDev.set(30.0);
4165 292 : doPBCorr = true;
4166 292 : conjBeams = true;
4167 292 : computePAStep=360.0;
4168 292 : rotatePAStep=5.0;
4169 :
4170 : // extra params for single-dish
4171 292 : pointingDirCol = "";
4172 292 : convertFirst = "never";
4173 292 : skyPosThreshold = 0.0;
4174 292 : convSupport = -1;
4175 292 : truncateSize = Quantity(-1.0);
4176 292 : gwidth = Quantity(-1.0);
4177 292 : jwidth = Quantity(-1.0);
4178 292 : minWeight = 0.0;
4179 292 : clipMinMax = False;
4180 :
4181 : // Mapper type
4182 292 : mType = String("default");
4183 :
4184 292 : }
4185 :
4186 0 : Record SynthesisParamsGrid::toRecord() const
4187 : {
4188 0 : Record gridpar;
4189 :
4190 0 : gridpar.define("imagename", imageName);
4191 : // FTMachine params
4192 0 : gridpar.define("padding", padding);
4193 0 : gridpar.define("useautocorr",useAutoCorr );
4194 0 : gridpar.define("usedoubleprec", useDoublePrec);
4195 0 : gridpar.define("wprojplanes", wprojplanes);
4196 0 : gridpar.define("convfunc", convFunc);
4197 0 : gridpar.define("vptable", vpTable);
4198 :
4199 0 : gridpar.define("facets", facets);
4200 0 : gridpar.define("chanchunks", chanchunks);
4201 :
4202 0 : gridpar.define("interpolation",interpolation);
4203 :
4204 0 : gridpar.define("distance", QuantityToString(distance));
4205 0 : gridpar.define("tracksource", trackSource);
4206 0 : gridpar.define("trackdir", MDirectionToString( trackDir ));
4207 :
4208 0 : gridpar.define("aterm",aTermOn );
4209 0 : gridpar.define("psterm",psTermOn );
4210 0 : gridpar.define("mterm",mTermOn );
4211 0 : gridpar.define("wbawp", wbAWP);
4212 0 : gridpar.define("cfcache", cfCache);
4213 0 : gridpar.define("usepointing",usePointing );
4214 0 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
4215 0 : gridpar.define("dopbcorr", doPBCorr);
4216 0 : gridpar.define("conjbeams",conjBeams );
4217 0 : gridpar.define("computepastep", computePAStep);
4218 0 : gridpar.define("rotatepastep", rotatePAStep);
4219 :
4220 0 : gridpar.define("pointingcolumntouse", pointingDirCol );
4221 0 : gridpar.define("convertfirst", convertFirst );
4222 0 : gridpar.define("skyposthreshold", skyPosThreshold );
4223 0 : gridpar.define("convsupport", convSupport );
4224 0 : gridpar.define("truncate", QuantityToString(truncateSize) );
4225 0 : gridpar.define("gwidth", QuantityToString(gwidth) );
4226 0 : gridpar.define("jwidth", QuantityToString(jwidth) );
4227 0 : gridpar.define("minweight", minWeight );
4228 0 : gridpar.define("clipminmax", clipMinMax );
4229 :
4230 0 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
4231 : /// else gridpar.define("deconvolver","singleterm");
4232 :
4233 0 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
4234 0 : else gridpar.define("gridder", gridder);
4235 :
4236 : // gridpar.define("mtype", mType);
4237 :
4238 0 : return gridpar;
4239 0 : }
4240 :
4241 :
4242 :
4243 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
4244 :
4245 : /////////////////////// Deconvolver Parameters
4246 :
4247 146 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4248 : {
4249 146 : setDefaults();
4250 146 : }
4251 :
4252 146 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4253 : {
4254 146 : }
4255 :
4256 :
4257 146 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4258 : {
4259 146 : setDefaults();
4260 :
4261 146 : String err("");
4262 :
4263 : try
4264 : {
4265 :
4266 146 : err += readVal( inrec, String("imagename"), imageName );
4267 146 : err += readVal( inrec, String("deconvolver"), algorithm );
4268 :
4269 :
4270 : //err += readVal( inrec, String("startmodel"), startModel );
4271 : // startmodel parsing copied from SynthesisParamsImage. Clean this up !!!
4272 146 : if( inrec.isDefined("startmodel") )
4273 : {
4274 146 : if( inrec.dataType("startmodel")==TpString )
4275 : {
4276 73 : String onemodel;
4277 73 : err += readVal( inrec, String("startmodel"), onemodel );
4278 73 : if( onemodel.length()>0 )
4279 : {
4280 0 : startModel.resize(1);
4281 0 : startModel[0] = onemodel;
4282 : }
4283 73 : else {startModel.resize();}
4284 73 : }
4285 146 : else if( inrec.dataType("startmodel")==TpArrayString ||
4286 73 : inrec.dataType("startmodel")==TpArrayBool)
4287 : {
4288 73 : err += readVal( inrec, String("startmodel"), startModel );
4289 : }
4290 : else {
4291 0 : err += String("startmodel must be either a string(singleterm) or a list of strings(multiterm)\n");
4292 : }
4293 : }
4294 : //------------------------
4295 :
4296 146 : err += readVal( inrec, String("id"), deconvolverId );
4297 146 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4298 :
4299 146 : err += readVal( inrec, String("scales"), scales );
4300 146 : err += readVal( inrec, String("scalebias"), scalebias );
4301 :
4302 146 : err += readVal( inrec, String("usemask"), maskType );
4303 146 : if( maskType=="auto-thresh" )
4304 : {
4305 0 : autoMaskAlgorithm = "thresh";
4306 : }
4307 146 : else if( maskType=="auto-thesh2" )
4308 : {
4309 0 : autoMaskAlgorithm = "thresh2";
4310 : }
4311 146 : else if( maskType=="auto-onebox" )
4312 : {
4313 0 : autoMaskAlgorithm = "onebox";
4314 : }
4315 146 : else if( maskType=="user" || maskType=="pb" )
4316 : {
4317 146 : autoMaskAlgorithm = "";
4318 : }
4319 :
4320 :
4321 146 : if( inrec.isDefined("mask") )
4322 : {
4323 146 : if( inrec.dataType("mask")==TpString )
4324 : {
4325 146 : err+= readVal( inrec, String("mask"), maskString );
4326 : }
4327 0 : else if( inrec.dataType("mask")==TpArrayString )
4328 : {
4329 0 : err+= readVal( inrec, String("mask"), maskList );
4330 : }
4331 : }
4332 :
4333 146 : if( inrec.isDefined("pbmask") )
4334 : {
4335 146 : err += readVal( inrec, String("pbmask"), pbMask );
4336 : }
4337 146 : if( inrec.isDefined("maskthreshold") )
4338 : {
4339 146 : if( inrec.dataType("maskthreshold")==TpString )
4340 : {
4341 146 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4342 : //deal with the case a string is a float value without unit
4343 146 : Quantity testThresholdString;
4344 146 : Quantity::read(testThresholdString,maskThreshold);
4345 146 : if( testThresholdString.getUnit()=="" )
4346 : {
4347 146 : if(testThresholdString.getValue()<1.0)
4348 : {
4349 146 : fracOfPeak = testThresholdString.getValue();
4350 146 : maskThreshold=String("");
4351 : }
4352 : }
4353 146 : }
4354 0 : else if( inrec.dataType("maskthreshold")==TpFloat || inrec.dataType("maskthreshold")==TpDouble )
4355 : {
4356 :
4357 0 : err += readVal( inrec, String("maskthreshold"), fracOfPeak );
4358 0 : if( fracOfPeak >=1.0 )
4359 : {
4360 : // maskthreshold is sigma ( * rms = threshold)
4361 : //
4362 0 : maskThreshold=String::toString(fracOfPeak);
4363 0 : fracOfPeak=0.0;
4364 : }
4365 : }
4366 : else
4367 : {
4368 0 : err += "maskthreshold must be a string, float, or double\n";
4369 : }
4370 : }
4371 146 : if( inrec.isDefined("maskresolution") )
4372 : {
4373 146 : if( inrec.dataType("maskresolution")==TpString )
4374 : {
4375 146 : err += readVal(inrec, String("maskresolution"), maskResolution );
4376 : //deal with the case a string is a float value without unit
4377 146 : Quantity testResolutionString;
4378 146 : Quantity::read(testResolutionString,maskResolution);
4379 146 : if( testResolutionString.getUnit()=="" )
4380 : {
4381 146 : maskResByBeam = testResolutionString.getValue();
4382 146 : maskResolution=String("");
4383 : }
4384 146 : }
4385 0 : else if( inrec.dataType("maskresolution")==TpFloat || inrec.dataType("maskresolution")==TpDouble )
4386 : {
4387 :
4388 0 : err += readVal( inrec, String("maskresolution"), maskResByBeam );
4389 : }
4390 : else
4391 : {
4392 0 : err += "maskresolution must be a string, float, or double\n";
4393 : }
4394 : }
4395 :
4396 :
4397 146 : if( inrec.isDefined("nmask") )
4398 : {
4399 146 : if( inrec.dataType("nmask")==TpInt )
4400 : {
4401 146 : err+= readVal(inrec, String("nmask"), nMask );
4402 : }
4403 : else
4404 : {
4405 0 : err+= "nmask must be an integer\n";
4406 : }
4407 : }
4408 146 : if( inrec.isDefined("autoadjust") )
4409 : {
4410 73 : if( inrec.dataType("autoadjust")==TpBool )
4411 : {
4412 73 : err+= readVal(inrec, String("autoadjust"), autoAdjust );
4413 : }
4414 : else
4415 : {
4416 0 : err+= "autoadjust must be a bool\n";
4417 : }
4418 : }
4419 : //param for the Asp-Clean to trigger Hogbom Clean
4420 146 : if (inrec.isDefined("fusedthreshold"))
4421 : {
4422 146 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4423 146 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4424 : else
4425 0 : err += "fusedthreshold must be a float or double";
4426 : }
4427 146 : if (inrec.isDefined("specmode"))
4428 : {
4429 146 : if(inrec.dataType("specmode") == TpString)
4430 146 : err += readVal(inrec, String("specmode"), specmode);
4431 : else
4432 0 : err += "specmode must be a string";
4433 : }
4434 : //largest scale size for the Asp-Clean to overwrite the default
4435 146 : if (inrec.isDefined("largestscale"))
4436 : {
4437 146 : if (inrec.dataType("largestscale") == TpInt)
4438 146 : err += readVal(inrec, String("largestscale"), largestscale);
4439 : else
4440 0 : err += "largestscale must be an integer";
4441 : }
4442 : //params for the new automasking algorithm
4443 146 : if( inrec.isDefined("sidelobethreshold"))
4444 : {
4445 146 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4446 : {
4447 146 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4448 : }
4449 : else
4450 : {
4451 0 : err+= "sidelobethreshold must be a float or double";
4452 : }
4453 : }
4454 :
4455 146 : if( inrec.isDefined("noisethreshold"))
4456 : {
4457 146 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4458 : {
4459 146 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4460 : }
4461 : else
4462 : {
4463 0 : err+= "noisethreshold must be a float or double";
4464 : }
4465 : }
4466 146 : if( inrec.isDefined("lownoisethreshold"))
4467 : {
4468 146 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4469 : {
4470 146 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4471 : }
4472 : else
4473 : {
4474 0 : err+= "lownoisethreshold must be a float or double";
4475 : }
4476 : }
4477 146 : if( inrec.isDefined("negativethreshold"))
4478 : {
4479 146 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4480 : {
4481 146 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4482 : }
4483 : else
4484 : {
4485 0 : err+= "negativethreshold must be a float or double";
4486 : }
4487 : }
4488 146 : if( inrec.isDefined("smoothfactor"))
4489 : {
4490 146 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4491 : {
4492 146 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4493 : }
4494 : else
4495 : {
4496 0 : err+= "smoothfactor must be a float or double";
4497 : }
4498 : }
4499 146 : if( inrec.isDefined("minbeamfrac"))
4500 : {
4501 146 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4502 : {
4503 146 : err+= readVal(inrec, String("minbeamfrac"), minBeamFrac );
4504 : }
4505 : else
4506 : {
4507 0 : if (inrec.dataType("minbeamfrac")==TpInt) {
4508 0 : cerr<<"minbeamfrac is int"<<endl;
4509 : }
4510 0 : if (inrec.dataType("minbeamfrac")==TpString) {
4511 0 : cerr<<"minbeamfrac is String"<<endl;
4512 : }
4513 0 : err+= "minbeamfrac must be a float or double";
4514 : }
4515 : }
4516 146 : if( inrec.isDefined("cutthreshold"))
4517 : {
4518 146 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4519 : {
4520 146 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4521 : }
4522 : else {
4523 0 : err+= "cutthreshold must be a float or double";
4524 : }
4525 : }
4526 146 : if( inrec.isDefined("growiterations"))
4527 : {
4528 146 : if (inrec.dataType("growiterations")==TpInt) {
4529 146 : err+= readVal(inrec, String("growiterations"), growIterations );
4530 : }
4531 : else {
4532 0 : err+= "growiterations must be an integer\n";
4533 : }
4534 : }
4535 146 : if( inrec.isDefined("dogrowprune"))
4536 : {
4537 146 : if (inrec.dataType("dogrowprune")==TpBool) {
4538 146 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4539 : }
4540 : else {
4541 0 : err+= "dogrowprune must be a bool\n";
4542 : }
4543 : }
4544 146 : if( inrec.isDefined("minpercentchange"))
4545 : {
4546 146 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4547 146 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4548 : }
4549 : else {
4550 0 : err+= "minpercentchange must be a float or double";
4551 : }
4552 : }
4553 146 : if( inrec.isDefined("verbose"))
4554 : {
4555 146 : if (inrec.dataType("verbose")==TpBool ) {
4556 146 : err+= readVal(inrec, String("verbose"), verbose);
4557 : }
4558 : else {
4559 0 : err+= "verbose must be a bool";
4560 : }
4561 : }
4562 146 : if( inrec.isDefined("fastnoise"))
4563 : {
4564 146 : if (inrec.dataType("fastnoise")==TpBool ) {
4565 146 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4566 : }
4567 : else {
4568 0 : err+= "fastnoise must be a bool";
4569 : }
4570 : }
4571 146 : if( inrec.isDefined("nsigma") )
4572 : {
4573 146 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4574 146 : err+= readVal(inrec, String("nsigma"), nsigma );
4575 : }
4576 0 : else if(inrec.dataType("nsigma")==TpInt)
4577 : {
4578 : int tnsigma;
4579 0 : err+= readVal(inrec, String("nsigma"), tnsigma );
4580 0 : nsigma = float(tnsigma);
4581 : }
4582 : else {
4583 0 : err+= "nsigma must be an int, float or double";
4584 : }
4585 : }
4586 146 : if( inrec.isDefined("noRequireSumwt") )
4587 : {
4588 73 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4589 73 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4590 : }
4591 : else {
4592 0 : err+= "noRequireSumwt must be a bool";
4593 : }
4594 : }
4595 146 : if( inrec.isDefined("fullsummary") )
4596 : {
4597 146 : if (inrec.dataType("fullsummary")==TpBool) {
4598 146 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4599 : }
4600 : else {
4601 0 : err+= "fullsummary must be a bool";
4602 : }
4603 : }
4604 146 : if( inrec.isDefined("restoringbeam") )
4605 : {
4606 73 : String errinfo("");
4607 : try {
4608 :
4609 73 : if( inrec.dataType("restoringbeam")==TpString )
4610 : {
4611 0 : err += readVal( inrec, String("restoringbeam"), usebeam);
4612 : // FIXME ! usebeam.length() == 0 is a poorly formed conditional, it
4613 : // probably needs simplification or parenthesis, the compiler is
4614 : // compaining about it
4615 0 : if( (! usebeam.matches("common")) && usebeam.length()!=0 )
4616 : {
4617 0 : Quantity bsize;
4618 0 : err += readVal( inrec, String("restoringbeam"), bsize );
4619 0 : restoringbeam.setMajorMinor( bsize, bsize );
4620 0 : usebeam = String("");
4621 0 : }
4622 0 : errinfo = usebeam;
4623 : }
4624 73 : else if( inrec.dataType("restoringbeam")==TpArrayString )
4625 : {
4626 0 : Vector<String> bpars;
4627 0 : err += readVal( inrec, String("restoringbeam"), bpars );
4628 :
4629 0 : for (uInt i=0;i<bpars.nelements();i++) { errinfo += bpars[i] + " "; }
4630 :
4631 0 : if( bpars.nelements()==1 && bpars[0].length()>0 ) {
4632 0 : if( bpars[0]=="common") { usebeam="common"; }
4633 : else {
4634 0 : Quantity axis; stringToQuantity( bpars[0] , axis);
4635 0 : restoringbeam.setMajorMinor( axis, axis );
4636 0 : }
4637 0 : }else if( bpars.nelements()==2 ) {
4638 0 : Quantity majaxis, minaxis;
4639 0 : stringToQuantity( bpars[0], majaxis ); stringToQuantity( bpars[1], minaxis );
4640 0 : restoringbeam.setMajorMinor( majaxis, minaxis );
4641 0 : }else if( bpars.nelements()==3 ) {
4642 0 : Quantity majaxis, minaxis, pa;
4643 0 : stringToQuantity( bpars[0], majaxis ); stringToQuantity( bpars[1], minaxis ); stringToQuantity( bpars[2], pa );
4644 0 : restoringbeam.setMajorMinor( majaxis, minaxis );
4645 0 : restoringbeam.setPA( pa );
4646 0 : }else {
4647 0 : restoringbeam = GaussianBeam();
4648 0 : usebeam = String("");
4649 : }
4650 0 : }
4651 0 : } catch( AipsError &x) {
4652 0 : err += "Cannot construct a restoringbeam from supplied parameters " + errinfo + ". Please check that majoraxis >= minoraxis and all entries are strings.";
4653 0 : restoringbeam = GaussianBeam();
4654 0 : usebeam = String("");
4655 0 : }
4656 :
4657 73 : }// if isdefined(restoringbeam)
4658 :
4659 146 : if( inrec.isDefined("interactive") )
4660 : {
4661 146 : if( inrec.dataType("interactive")==TpBool )
4662 146 : {err += readVal( inrec, String("interactive"), interactive );}
4663 0 : else if ( inrec.dataType("interactive")==TpInt )
4664 0 : {Int inter=0; err += readVal( inrec, String("interactive"), inter); interactive=(Bool)inter;}
4665 : }
4666 :
4667 : //err += readVal( inrec, String("interactive"), interactive );
4668 :
4669 146 : err += verify();
4670 :
4671 : }
4672 0 : catch(AipsError &x)
4673 : {
4674 0 : err = err + x.getMesg() + "\n";
4675 0 : }
4676 :
4677 146 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4678 :
4679 146 : }
4680 :
4681 146 : String SynthesisParamsDeconv::verify() const
4682 : {
4683 146 : String err;
4684 :
4685 146 : if( imageName=="" ) {err += "Please supply an image name\n";}
4686 :
4687 : // Allow mask inputs in only one way. User specified OR already on disk. Not both
4688 146 : if( maskString.length()>0 )
4689 : {
4690 0 : File fp( imageName+".mask" );
4691 0 : if( fp.exists() ) err += "Mask image " + imageName+".mask exists, but a specific input mask of " + maskString + " has also been supplied. Please either reset mask='' to reuse the existing mask, or delete " + imageName + ".mask before restarting";
4692 0 : }
4693 :
4694 146 : if( pbMask >= 1.0)
4695 0 : {err += "pbmask must be < 1.0 \n"; }
4696 146 : else if( pbMask < 0.0)
4697 0 : {err += "pbmask must be a positive value \n"; }
4698 :
4699 146 : if( maskType=="none" )
4700 : {
4701 0 : if( maskString!="" || (maskList.nelements()!=0 && maskList[0]!="") )
4702 : {
4703 0 : cerr<<"maskString="<<maskString<<endl;
4704 0 : cerr<<"maskList.nelements()="<<maskList.nelements()<<" maskList[0]="<<maskList[0]<<endl;
4705 0 : err += "mask is specified but usemask='none'. Please set usemask='user' to use the mask parameter\n";}
4706 : }
4707 146 : if ( fracOfPeak >= 1.0)
4708 0 : {err += "fracofpeak must be < 1.0 \n"; }
4709 146 : else if ( fracOfPeak < 0.0)
4710 0 : {err += "fracofpeak must be a positive value \n"; }
4711 :
4712 146 : return err;
4713 0 : }
4714 :
4715 292 : void SynthesisParamsDeconv::setDefaults()
4716 : {
4717 292 : imageName="";
4718 292 : algorithm="hogbom";
4719 292 : startModel=Vector<String>(0);
4720 292 : deconvolverId=0;
4721 292 : nTaylorTerms=1;
4722 292 : scales.resize(1); scales[0]=0.0;
4723 292 : scalebias=0.6;
4724 292 : maskType="none";
4725 292 : maskString="";
4726 292 : maskList.resize(1); maskList[0]="";
4727 292 : pbMask=0.0;
4728 292 : autoMaskAlgorithm="thresh";
4729 292 : maskThreshold="";
4730 292 : maskResolution="";
4731 292 : fracOfPeak=0.0;
4732 292 : nMask=0;
4733 292 : interactive=false;
4734 292 : autoAdjust=False;
4735 292 : fusedThreshold = 0.0;
4736 292 : specmode="mfs";
4737 292 : largestscale = -1;
4738 292 : }
4739 :
4740 73 : Record SynthesisParamsDeconv::toRecord() const
4741 : {
4742 73 : Record decpar;
4743 :
4744 73 : decpar.define("imagename", imageName);
4745 73 : decpar.define("deconvolver", algorithm);
4746 73 : decpar.define("startmodel",startModel);
4747 73 : decpar.define("id",deconvolverId);
4748 73 : decpar.define("nterms",nTaylorTerms);
4749 73 : decpar.define("scales",scales);
4750 73 : decpar.define("scalebias",scalebias);
4751 73 : decpar.define("usemask",maskType);
4752 73 : decpar.define("fusedthreshold", fusedThreshold);
4753 73 : decpar.define("specmode", specmode);
4754 73 : decpar.define("largestscale", largestscale);
4755 73 : if( maskList.nelements()==1 && maskList[0]=="")
4756 : {
4757 73 : decpar.define("mask",maskString);
4758 : }
4759 : else {
4760 0 : decpar.define("mask",maskList);
4761 : }
4762 73 : decpar.define("pbmask",pbMask);
4763 73 : if (fracOfPeak > 0.0)
4764 : {
4765 0 : decpar.define("maskthreshold",fracOfPeak);
4766 : }
4767 : else
4768 : {
4769 73 : decpar.define("maskthreshold",maskThreshold);
4770 : }
4771 73 : decpar.define("maskresolution",maskResolution);
4772 73 : decpar.define("nmask",nMask);
4773 73 : decpar.define("autoadjust",autoAdjust);
4774 73 : decpar.define("sidelobethreshold",sidelobeThreshold);
4775 73 : decpar.define("noisethreshold",noiseThreshold);
4776 73 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4777 73 : decpar.define("negativethreshold",negativeThreshold);
4778 73 : decpar.define("smoothfactor",smoothFactor);
4779 73 : decpar.define("minbeamfrac",minBeamFrac);
4780 73 : decpar.define("cutthreshold",cutThreshold);
4781 73 : decpar.define("growiterations",growIterations);
4782 73 : decpar.define("dogrowprune",doGrowPrune);
4783 73 : decpar.define("minpercentchange",minPercentChange);
4784 73 : decpar.define("verbose", verbose);
4785 73 : decpar.define("fastnoise", fastnoise);
4786 73 : decpar.define("interactive",interactive);
4787 73 : decpar.define("nsigma",nsigma);
4788 73 : decpar.define("noRequireSumwt",noRequireSumwt);
4789 73 : decpar.define("fullsummary",fullsummary);
4790 :
4791 73 : return decpar;
4792 0 : }
4793 :
4794 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4795 :
4796 :
4797 : } //# NAMESPACE CASA - END
4798 :
|
