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 2 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 2 : }
89 :
90 2 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 2 : }
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 0 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 0 : Int N=vb.nRows(),M=-1;
108 0 : for(Int i=0;i<N;i++)
109 : {
110 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 0 : {M++;break;}
112 : }
113 0 : 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 0 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 0 : Int n=npix;
121 :
122 0 : if( n%2 !=0 ){ n+= 1; }
123 :
124 0 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 0 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 0 : 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 0 : newlarge=product(fac);
136 0 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 0 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 0 : return newlarge;
141 0 : }
142 :
143 : // Return the prime factors of the given number
144 0 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 0 : Vector<uInt> factors;
147 :
148 0 : Int lastresult = n;
149 0 : Int sqlast = int(sqrt(n))+1;
150 :
151 0 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 0 : Int c=2;
153 : while(1)
154 : {
155 0 : if( lastresult == 1 || c > sqlast ) { break; }
156 0 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 0 : if( c>sqlast){ c=lastresult; break; }
160 0 : if( lastresult % c == 0 ) { break; }
161 0 : c += 1;
162 : }
163 0 : factors.resize( factors.nelements()+1, true );
164 0 : factors[ factors.nelements()-1 ] = c;
165 0 : lastresult /= c;
166 : }
167 0 : 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 0 : 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 0 : 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 0 : Record outRec;
705 0 : outRec.define("type", type);
706 0 : outRec.define("rmode", rmode);
707 0 : Record quantRec;
708 0 : QuantumHolder(noise).toRecord(quantRec);
709 0 : outRec.defineRecord("noise", quantRec);
710 0 : outRec.define("robust", robust);
711 0 : QuantumHolder(fieldofview).toRecord(quantRec);
712 0 : outRec.defineRecord("fieldofview", quantRec);
713 0 : outRec.define("npixels", npixels);
714 0 : outRec.define("multifield", multiField);
715 0 : outRec.define("usecubebriggs", useCubeBriggs);
716 0 : outRec.define("filtertype", filtertype);
717 0 : QuantumHolder(filterbmaj).toRecord(quantRec);
718 0 : outRec.defineRecord("filterbmaj", quantRec);
719 0 : QuantumHolder(filterbmin).toRecord(quantRec);
720 0 : outRec.defineRecord("filterbmin", quantRec);
721 0 : QuantumHolder(filterbpa).toRecord(quantRec);
722 0 : outRec.defineRecord("filterbpa", quantRec);
723 0 : outRec.define("fracBW", fracBW);
724 :
725 0 : return outRec;
726 0 : }
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 0 : 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 0 : Bool isOn = false;
924 0 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
925 0 : if (!isOn)
926 0 : 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 0 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1421 : {
1422 0 : if( rec.isDefined( id ) )
1423 : {
1424 0 : String inval("");
1425 0 : if( rec.dataType( id )==TpString )
1426 0 : { 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 0 : if(inval.length()>0){val=inval;}
1432 0 : return String("");
1433 : }
1434 0 : else { return String(id + " must be a string\n"); }
1435 0 : }
1436 0 : else { return String("");}
1437 : }
1438 :
1439 : // Read Integer from Record
1440 0 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1441 : {
1442 0 : if( rec.isDefined( id ) )
1443 : {
1444 0 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1445 0 : else { return String(id + " must be an integer\n"); }
1446 : }
1447 0 : else { return String(""); }
1448 : }
1449 :
1450 : // Read Float from Record
1451 0 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1452 : {
1453 0 : if( rec.isDefined( id ) )
1454 : {
1455 0 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1456 0 : { rec.get( id , val ); return String(""); }
1457 0 : else { return String(id + " must be a float\n"); }
1458 : }
1459 0 : else { return String(""); }
1460 : }
1461 :
1462 : // Read Bool from Record
1463 0 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1464 : {
1465 0 : if( rec.isDefined( id ) )
1466 : {
1467 0 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1468 0 : else { return String(id + " must be a bool\n"); }
1469 : }
1470 0 : 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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1493 : {
1494 0 : if( rec.isDefined( id ) )
1495 : {
1496 0 : if( rec.dataType( id )==TpArrayFloat )
1497 : {
1498 0 : 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 0 : 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 0 : 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 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1526 : {
1527 0 : Vector<Bool> vec; rec.get( id, vec );
1528 0 : 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 0 : }
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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1539 : {
1540 0 : if( rec.isDefined( id ) )
1541 : {
1542 0 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1543 0 : { 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 0 : 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 0 : 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 0 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1570 : {
1571 : try
1572 : {
1573 0 : String tmpRF, tmpRA, tmpDEC;
1574 0 : std::istringstream iss(instr);
1575 0 : iss >> tmpRF >> tmpRA >> tmpDEC;
1576 0 : 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 0 : casacore::Quantity tmpQRA;
1583 0 : casacore::Quantity tmpQDEC;
1584 0 : casacore::Quantity::read(tmpQRA, tmpRA);
1585 0 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1586 :
1587 0 : 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 0 : Bool status = MDirection::getType(theRF, tmpRF);
1598 0 : if (!status) {
1599 0 : throw AipsError();
1600 : }
1601 0 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1602 0 : return String("");
1603 0 : }
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 0 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1612 : {
1613 0 : if( rec.isDefined( id ) )
1614 : {
1615 0 : if( rec.dataType( id )==TpString )
1616 0 : { 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 0 : else{ return String(""); }
1620 : }
1621 :
1622 : // Read MDirection from a Record string
1623 0 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1624 : {
1625 0 : if( rec.isDefined( id ) )
1626 : {
1627 0 : if( rec.dataType( id )==TpString )
1628 0 : { 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 0 : 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 0 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1683 : {
1684 0 : setDefaults();
1685 0 : }
1686 :
1687 0 : SynthesisParamsSelect::~SynthesisParamsSelect()
1688 : {
1689 0 : }
1690 :
1691 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1692 0 : operator=(other);
1693 0 : }
1694 :
1695 0 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1696 0 : if(this!=&other) {
1697 0 : msname=other.msname;
1698 0 : spw=other.spw;
1699 0 : freqbeg=other.freqbeg;
1700 0 : freqend=other.freqend;
1701 0 : freqframe=other.freqframe;
1702 0 : field=other.field;
1703 0 : antenna=other.antenna;
1704 0 : timestr=other.timestr;
1705 0 : scan=other.scan;
1706 0 : obs=other.obs;
1707 0 : state=other.state;
1708 0 : uvdist=other.uvdist;
1709 0 : taql=other.taql;
1710 0 : usescratch=other.usescratch;
1711 0 : readonly=other.readonly;
1712 0 : incrmodel=other.incrmodel;
1713 0 : datacolumn=other.datacolumn;
1714 :
1715 : }
1716 0 : return *this;
1717 : }
1718 :
1719 0 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1720 : {
1721 0 : setDefaults();
1722 :
1723 0 : String err("");
1724 :
1725 : try
1726 : {
1727 :
1728 0 : err += readVal( inrec, String("msname"), msname );
1729 :
1730 0 : err += readVal( inrec, String("readonly"), readonly );
1731 0 : err += readVal( inrec, String("usescratch"), usescratch );
1732 :
1733 : // override with entries from savemodel.
1734 0 : String savemodel("");
1735 0 : err += readVal( inrec, String("savemodel"), savemodel );
1736 0 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1737 0 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1738 0 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1739 :
1740 0 : err += readVal( inrec, String("incrmodel"), incrmodel );
1741 :
1742 0 : err += readVal( inrec, String("spw"), spw );
1743 :
1744 : /// -------------------------------------------------------------------------------------
1745 : // Why are these params here ? Repeats of defineimage.
1746 0 : err += readVal( inrec, String("freqbeg"), freqbeg);
1747 0 : err += readVal( inrec, String("freqend"), freqend);
1748 :
1749 0 : String freqframestr( MFrequency::showType(freqframe) );
1750 0 : err += readVal( inrec, String("outframe"), freqframestr);
1751 0 : if( ! MFrequency::getType(freqframe, freqframestr) )
1752 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1753 : /// -------------------------------------------------------------------------------------
1754 :
1755 0 : err += readVal( inrec, String("field"),field);
1756 0 : err += readVal( inrec, String("antenna"),antenna);
1757 0 : err += readVal( inrec, String("timestr"),timestr);
1758 0 : err += readVal( inrec, String("scan"),scan);
1759 0 : err += readVal( inrec, String("obs"),obs);
1760 0 : err += readVal( inrec, String("state"),state);
1761 0 : err += readVal( inrec, String("uvdist"),uvdist);
1762 0 : err += readVal( inrec, String("taql"),taql);
1763 :
1764 0 : err += readVal( inrec, String("datacolumn"),datacolumn);
1765 :
1766 0 : err += verify();
1767 :
1768 0 : }
1769 0 : catch(AipsError &x)
1770 : {
1771 0 : err = err + x.getMesg() + "\n";
1772 0 : }
1773 :
1774 0 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1775 :
1776 0 : }
1777 :
1778 0 : String SynthesisParamsSelect::verify() const
1779 : {
1780 0 : String err;
1781 : // Does the MS exist on disk.
1782 0 : Directory thems( msname );
1783 0 : if( thems.exists() )
1784 : {
1785 : // Is it readable ?
1786 0 : if( ! thems.isReadable() )
1787 0 : { err += "MS " + msname + " is not readable.\n"; }
1788 : // Depending on 'readonly', is the MS writable ?
1789 0 : 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 0 : return err;
1796 0 : }
1797 :
1798 0 : void SynthesisParamsSelect::setDefaults()
1799 : {
1800 0 : msname="";
1801 0 : spw="";
1802 0 : freqbeg="";
1803 0 : freqend="";
1804 0 : MFrequency::getType(freqframe,"LSRK");
1805 0 : field="";
1806 0 : antenna="";
1807 0 : timestr="";
1808 0 : scan="";
1809 0 : obs="";
1810 0 : state="";
1811 0 : uvdist="";
1812 0 : taql="";
1813 0 : usescratch=false;
1814 0 : readonly=true;
1815 0 : incrmodel=false;
1816 0 : datacolumn="corrected";
1817 0 : }
1818 :
1819 0 : Record SynthesisParamsSelect::toRecord()const
1820 : {
1821 0 : Record selpar;
1822 0 : selpar.define("msname",msname);
1823 0 : selpar.define("spw",spw);
1824 0 : selpar.define("freqbeg",freqbeg);
1825 0 : selpar.define("freqend",freqend);
1826 0 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1827 : //looks like fromRecord looks for outframe !
1828 0 : selpar.define("outframe", MFrequency::showType(freqframe));
1829 0 : selpar.define("field",field);
1830 0 : selpar.define("antenna",antenna);
1831 0 : selpar.define("timestr",timestr);
1832 0 : selpar.define("scan",scan);
1833 0 : selpar.define("obs",obs);
1834 0 : selpar.define("state",state);
1835 0 : selpar.define("uvdist",uvdist);
1836 0 : selpar.define("taql",taql);
1837 0 : selpar.define("usescratch",usescratch);
1838 0 : selpar.define("readonly",readonly);
1839 0 : selpar.define("incrmodel",incrmodel);
1840 0 : selpar.define("datacolumn",datacolumn);
1841 :
1842 0 : return selpar;
1843 0 : }
1844 :
1845 :
1846 : /////////////////////// Image Parameters
1847 :
1848 0 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1849 : {
1850 0 : setDefaults();
1851 0 : }
1852 :
1853 0 : SynthesisParamsImage::~SynthesisParamsImage()
1854 : {
1855 0 : }
1856 :
1857 0 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1858 0 : if(this != &other){
1859 0 : imageName=other.imageName;
1860 0 : stokes=other.stokes;
1861 0 : startModel.resize(); startModel=other.startModel;
1862 0 : imsize.resize(); imsize=other.imsize;
1863 0 : cellsize.resize(); cellsize=other.cellsize;
1864 0 : projection=other.projection;
1865 0 : useNCP=other.useNCP;
1866 0 : phaseCenter=other.phaseCenter;
1867 0 : phaseCenterFieldId=other.phaseCenterFieldId;
1868 0 : obslocation=other.obslocation;
1869 0 : pseudoi=other.pseudoi;
1870 0 : nchan=other.nchan;
1871 0 : nTaylorTerms=other.nTaylorTerms;
1872 0 : chanStart=other.chanStart;
1873 0 : chanStep=other.chanStep;
1874 0 : freqStart=other.freqStart;
1875 0 : freqStep=other.freqStep;
1876 0 : refFreq=other.refFreq;
1877 0 : velStart=other.velStart;
1878 0 : velStep=other.velStep;
1879 0 : freqFrame=other.freqFrame;
1880 0 : mFreqStart=other.mFreqStart;
1881 0 : mFreqStep=other.mFreqStep;
1882 0 : mVelStart=other.mVelStart;
1883 0 : mVelStep=other.mVelStep;
1884 0 : restFreq.resize(); restFreq=other.restFreq;
1885 0 : start=other.start;
1886 0 : step=other.step;
1887 0 : frame=other.frame;
1888 0 : veltype=other.veltype;
1889 0 : mode=other.mode;
1890 0 : reffreq=other.reffreq;
1891 0 : sysvel=other.sysvel;
1892 0 : sysvelframe=other.sysvelframe;
1893 0 : sysvelvalue=other.sysvelvalue;
1894 0 : qmframe=other.qmframe;
1895 0 : mveltype=other.mveltype;
1896 0 : tststr=other.tststr;
1897 0 : startRecord=other.startRecord;
1898 0 : stepRecord=other.stepRecord;
1899 0 : reffreqRecord=other.reffreqRecord;
1900 0 : sysvelRecord=other.sysvelRecord;
1901 0 : restfreqRecord=other.restfreqRecord;
1902 0 : csysRecord=other.csysRecord;
1903 0 : csys=other.csys;
1904 0 : imshape.resize(); imshape=other.imshape;
1905 0 : freqFrameValid=other.freqFrameValid;
1906 0 : overwrite=other.overwrite;
1907 0 : deconvolver=other.deconvolver;
1908 0 : distance=other.distance;
1909 0 : trackDir=other.trackDir;
1910 0 : trackSource=other.trackSource;
1911 0 : movingSource=other.movingSource;
1912 :
1913 :
1914 :
1915 : }
1916 :
1917 0 : return *this;
1918 :
1919 : }
1920 :
1921 0 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec) {
1922 0 : fromRecord(inrec, False);
1923 0 : }
1924 :
1925 0 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec,
1926 : const casacore::Bool isSingleDish) {
1927 0 : setDefaults();
1928 0 : String err("");
1929 0 : LogIO os(LogOrigin("SynthesisParamsImage", "fromRecord"));
1930 :
1931 : try {
1932 0 : err += readVal(inrec, String("imagename"), imageName);
1933 :
1934 : //// imsize
1935 0 : if (inrec.isDefined("imsize")) {
1936 0 : DataType tp = inrec.dataType("imsize");
1937 0 : if (tp == TpInt) { // A single integer for both dimensions.
1938 : Int npix;
1939 0 : inrec.get("imsize", npix);
1940 0 : imsize.resize(2);
1941 0 : 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 0 : if (inrec.isDefined("cell")) {
1960 : try {
1961 0 : DataType tp = inrec.dataType("cell");
1962 0 : if (tp == TpInt ||
1963 0 : tp == TpFloat ||
1964 : tp == TpDouble) {
1965 0 : Double cell = inrec.asDouble("cell");
1966 0 : cellsize.set(Quantity(cell, "arcsec"));
1967 0 : } else if (tp == TpArrayInt ||
1968 0 : 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 0 : } else if (tp == TpString) {
1981 0 : String cell;
1982 0 : inrec.get("cell", cell);
1983 0 : Quantity qcell;
1984 0 : err += stringToQuantity(cell, qcell);
1985 0 : cellsize.set(qcell);
1986 0 : } 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 0 : err += readVal(inrec, String("stokes"), stokes);
2010 0 : if (stokes.matches("pseudoI")) {
2011 0 : stokes = "I";
2012 0 : pseudoi = true;
2013 : } else {
2014 0 : pseudoi = false;
2015 : }
2016 :
2017 : /// PseudoI
2018 :
2019 : ////nchan
2020 0 : err += readVal(inrec, String("nchan"), nchan);
2021 :
2022 : /// phaseCenter (as a string) . // Add INT support later.
2023 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2024 0 : if (inrec.isDefined("phasecenter")) {
2025 0 : String pcent("");
2026 0 : if (inrec.dataType("phasecenter") == TpString) {
2027 0 : inrec.get("phasecenter", pcent);
2028 0 : if (pcent.length() > 0) { // if it's zero length, it means 'figure out from first field in MS'.
2029 0 : err += readVal(inrec, String("phasecenter"), phaseCenter);
2030 0 : 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 0 : err += readVal(inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2034 : } else {
2035 0 : phaseCenterFieldId = -2; // deal with this later in buildCoordinateSystem to assign the first selected field
2036 : }
2037 : }
2038 0 : if (inrec.dataType("phasecenter") == TpInt ||
2039 0 : inrec.dataType("phasecenter") == TpFloat ||
2040 0 : 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 0 : if ((inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter") != TpInt
2045 0 : && inrec.dataType("phasecenter") != TpFloat && inrec.dataType("phasecenter") != TpDouble)) {
2046 0 : err += String("Cannot set phasecenter. Please specify a string or int\n");
2047 : }
2048 0 : }
2049 :
2050 : //// Projection
2051 0 : if (inrec.isDefined("projection")) {
2052 0 : if (inrec.dataType("projection") == TpString) {
2053 0 : String pstr;
2054 0 : inrec.get("projection", pstr);
2055 : try {
2056 0 : if (pstr.matches("NCP")) {
2057 0 : pstr = "SIN";
2058 0 : useNCP = true;
2059 : }
2060 0 : projection = Projection::type(pstr);
2061 0 : } catch (AipsError &x) {
2062 0 : err += String("Invalid projection code : " + pstr);
2063 0 : }
2064 0 : } else {
2065 0 : err += "projection must be a string\n";
2066 : }
2067 : } //projection
2068 :
2069 : // Frequency frame stuff.
2070 0 : err += readVal(inrec, String("specmode"), mode);
2071 : // Alias for 'mfs' is 'cont'
2072 0 : if (mode == "cont") mode = "mfs";
2073 :
2074 0 : err += readVal(inrec, String("outframe"), frame);
2075 0 : qmframe = "";
2076 : // mveltype is only set when start/step is given in mdoppler
2077 0 : mveltype = "";
2078 : //start
2079 0 : String startType("");
2080 0 : String widthType("");
2081 0 : if (inrec.isDefined("start")) {
2082 0 : if (inrec.dataType("start") == TpInt) {
2083 0 : err += readVal(inrec, String("start"), chanStart);
2084 0 : start = String::toString(chanStart);
2085 0 : startType = "chan";
2086 0 : } else if (inrec.dataType("start") == TpString) {
2087 0 : err += readVal(inrec, String("start"), start);
2088 0 : if (start.contains("Hz")) {
2089 0 : stringToQuantity(start, freqStart);
2090 0 : startType = "freq";
2091 0 : } 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 0 : if (inrec.isDefined("width")) {
2152 0 : if (inrec.dataType("width") == TpInt) {
2153 0 : err += readVal(inrec, String("width"), chanStep);
2154 0 : step = String::toString(chanStep);
2155 0 : widthType = "chan";
2156 0 : } else if (inrec.dataType("width") == TpString) {
2157 0 : err += readVal(inrec, String("width"), step);
2158 0 : if (step.contains("Hz")) {
2159 0 : stringToQuantity(step, freqStep);
2160 0 : widthType = "freq";
2161 0 : } 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 0 : 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 0 : if (inrec.isDefined("reffreq")) {
2226 0 : if (inrec.dataType("reffreq") == TpString) {
2227 0 : 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 0 : err += readVal(inrec, String("veltype"), veltype);
2254 0 : veltype = mveltype != "" ? mveltype : veltype;
2255 : // sysvel (String, Quantity)
2256 0 : if (inrec.isDefined("sysvel")) {
2257 0 : if (inrec.dataType("sysvel") == TpString) {
2258 0 : 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 0 : err += readVal(inrec, String("sysvelframe"), sysvelframe);
2270 :
2271 : // rest frequencies (record or vector of Strings)
2272 0 : if (inrec.isDefined("restfreq")) {
2273 0 : Vector<String> rfreqs(0);
2274 0 : Record restfreqSubRecord;
2275 0 : 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 0 : } else if (inrec.dataType("restfreq") == TpArrayString) {
2286 0 : err += readVal(inrec, String("restfreq"), rfreqs);
2287 0 : } else if (inrec.dataType("restfreq") == TpString) {
2288 0 : rfreqs.resize(1);
2289 0 : err += readVal(inrec, String("restfreq"), rfreqs[0]);
2290 : }
2291 0 : restFreq.resize(rfreqs.nelements());
2292 0 : for (uInt fr = 0; fr < rfreqs.nelements(); fr++) {
2293 0 : err += stringToQuantity(rfreqs[fr], restFreq[fr]);
2294 : }
2295 0 : } // if def restfreq
2296 :
2297 : // optional - coordsys, imshape
2298 : // if exist use them. May need a consistency check with the rest of impars?
2299 0 : 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 0 : String freqframestr = (frame != "" && qmframe != "") ? qmframe : frame;
2318 0 : if (frame != "" && !MFrequency::getType(freqFrame, freqframestr)) {
2319 0 : err += "Invalid Frequency Frame " + freqframestr;
2320 : }
2321 0 : err += readVal(inrec, String("restart"), overwrite);
2322 :
2323 0 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2324 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2325 0 : if (inrec.isDefined("startmodel")) {
2326 0 : if (inrec.dataType("startmodel") == TpString) {
2327 0 : String onemodel;
2328 0 : err += readVal(inrec, String("startmodel"), onemodel);
2329 0 : if (onemodel.length() > 0) {
2330 0 : startModel.resize(1);
2331 0 : startModel[0] = onemodel;
2332 : } else {
2333 0 : startModel.resize();
2334 : }
2335 0 : } 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 0 : err += readVal(inrec, String("nterms"), nTaylorTerms);
2344 0 : err += readVal(inrec, String("deconvolver"), deconvolver);
2345 :
2346 : // Force nchan=1 for anything other than cube modes...
2347 0 : if (mode == "mfs") nchan = 1;
2348 : //read obslocation
2349 0 : 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 0 : err += verify(isSingleDish);
2359 0 : }
2360 0 : catch (AipsError &x) {
2361 0 : err = err + x.getMesg() + "\n";
2362 0 : }
2363 :
2364 0 : if (err.length() > 0) throw(AipsError("Invalid Image Parameter set : " + err));
2365 0 : }
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 0 : String SynthesisParamsImage::verify(const casacore::Bool isSingleDish) const
2401 : {
2402 0 : LogIO os;
2403 0 : String err;
2404 :
2405 0 : if(imageName=="") {err += "Please supply an image name\n";}
2406 :
2407 0 : if(imsize.nelements() != 2){err += "imsize must be a vector of 2 Ints\n";}
2408 0 : if(cellsize.nelements() != 2) {err += "cellsize must be a vector of 2 Quantities\n";}
2409 0 : 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 0 : if((mode=="mfs") && nchan > 1){
2418 0 : err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";
2419 : }
2420 :
2421 0 : 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 0 : ! 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 0 : 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 0 : Int imxnew = SynthesisUtilMethods::getOptimumSize(imsize[0]);
2468 0 : Int imynew = SynthesisUtilMethods::getOptimumSize(imsize[1]);
2469 0 : 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 0 : return err;
2483 0 : }// 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 0 : void SynthesisParamsImage::setDefaults()
2496 : {
2497 : // Image definition parameters
2498 0 : imageName = String("");
2499 0 : imsize.resize(2); imsize.set(100);
2500 0 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2501 0 : stokes="I";
2502 0 : phaseCenter=MDirection();
2503 0 : phaseCenterFieldId=-1;
2504 0 : projection=Projection::SIN;
2505 0 : useNCP=false;
2506 0 : startModel=Vector<String>(0);
2507 0 : freqFrameValid=True;
2508 0 : overwrite=false;
2509 : // PseudoI
2510 0 : pseudoi=false;
2511 :
2512 : // Spectral coordinates
2513 0 : nchan=1;
2514 0 : mode="mfs";
2515 0 : start="";
2516 0 : step="";
2517 0 : chanStart=0;
2518 0 : 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 0 : freqStart=Quantity(0,"");
2524 0 : freqStep=Quantity(0,"");
2525 0 : velStart=Quantity(0,"");
2526 0 : velStep=Quantity(0,"");
2527 0 : veltype=String("radio");
2528 0 : restFreq.resize(0);
2529 0 : refFreq = Quantity(0,"Hz");
2530 0 : frame = "";
2531 0 : freqFrame=MFrequency::LSRK;
2532 0 : sysvel="";
2533 0 : sysvelframe="";
2534 0 : sysvelvalue=Quantity(0.0,"m/s");
2535 0 : nTaylorTerms=1;
2536 0 : deconvolver="hogbom";
2537 : ///csysRecord=Record();
2538 0 : }
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 0 : 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 0 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2674 0 : vi2.originChunks();
2675 0 : vi2.origin();
2676 : /// This version uses the new vi2/vb2
2677 : // get the first ms for multiple MSes
2678 : //MeasurementSet msobj=vi2.ms();
2679 0 : Int fld=vb->fieldId()(0);
2680 :
2681 : //handling first ms only
2682 0 : Double gfreqmax=-1.0;
2683 0 : Double gdatafend=-1.0;
2684 0 : Double gdatafstart=1e14;
2685 0 : Double gfreqmin=1e14;
2686 0 : Vector<Int> spwids0;
2687 0 : Int j=0;
2688 0 : Int minfmsid=0;
2689 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2690 0 : Double imStartFreq=getCubeImageStartFreq();
2691 0 : std::vector<Int> sourceMsWithStartFreq;
2692 :
2693 :
2694 0 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2695 : //auto forMS0=chansel.find(0);
2696 0 : map<Int, Vector<Int> > spwsels=forMS0->second;
2697 0 : Int nspws=spwsels.size();
2698 0 : Vector<Int> spwids(nspws);
2699 0 : Vector<Int> nChannels(nspws);
2700 0 : Vector<Int> firstChannels(nspws);
2701 : //Vector<Int> channelIncrement(nspws);
2702 :
2703 0 : Int k=0;
2704 0 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2705 0 : spwids[k]=it->first;
2706 0 : nChannels[k]=(it->second)[0];
2707 0 : firstChannels[k]=(it->second)[1];
2708 : }
2709 0 : if(j==0) {
2710 0 : spwids0.resize();
2711 0 : 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 0 : Double freqmin=0, freqmax=0;
2723 0 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2724 :
2725 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2726 0 : 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 0 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2732 : //cerr << "after " << datafstart << " " << datafend << endl;
2733 0 : if((datafstart > datafend) || !status)
2734 0 : throw(AipsError("spw selection failed"));
2735 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2736 :
2737 0 : if (mode=="cubedata") {
2738 0 : freqmin = datafstart;
2739 0 : freqmax = datafend;
2740 : }
2741 0 : 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 0 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2766 0 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2767 : //cerr << "after " << freqmin << " " << freqmax << endl;
2768 : }
2769 :
2770 :
2771 :
2772 :
2773 0 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2774 0 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2775 0 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2776 0 : 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 0 : 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 0 : }
2802 0 : 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 0 : MeasurementSet msobj = *mss[minfmsid];
2809 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2810 0 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2811 0 : }
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 0 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2844 : MeasurementSet& msobj,
2845 : Vector<Int> spwids, Int fld,
2846 : Double freqmin, Double freqmax,
2847 : Double datafstart, Double datafend )
2848 : {
2849 0 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2850 :
2851 0 : CoordinateSystem csys;
2852 0 : 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 0 : MSColumns msc(msobj);
2865 0 : String telescop = msc.observation().telescopeName()(0);
2866 0 : MEpoch obsEpoch = msc.timeMeas()(0);
2867 0 : MPosition obsPosition;
2868 0 : 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 0 : MDirection phaseCenterToUse = phaseCenter;
2880 :
2881 0 : if( phaseCenterFieldId != -1 )
2882 : {
2883 0 : MSFieldColumns msfield(msobj.field());
2884 0 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2885 : {
2886 0 : if(trackSource){
2887 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2888 : }
2889 : else{
2890 0 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2891 : }
2892 : }
2893 : else
2894 : {
2895 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2896 : }
2897 0 : }
2898 : // Setup Phase center if it is specified only by field id.
2899 :
2900 : /////////////////// Direction Coordinates
2901 0 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2902 : // Normalize correctly
2903 0 : MVAngle ra=mvPhaseCenter.get()(0);
2904 0 : ra(0.0);
2905 :
2906 0 : MVAngle dec=mvPhaseCenter.get()(1);
2907 0 : Vector<Double> refCoord(2);
2908 0 : refCoord(0)=ra.get().getValue();
2909 0 : refCoord(1)=dec;
2910 0 : Vector<Double> refPixel(2);
2911 0 : refPixel(0) = Double(imsize[0]/2);
2912 0 : refPixel(1) = Double(imsize[1]/2);
2913 :
2914 0 : Vector<Double> deltas(2);
2915 0 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2916 0 : deltas(1)=cellsize[1].get("rad").getValue();
2917 0 : Matrix<Double> xform(2,2);
2918 0 : xform=0.0;xform.diagonal()=1.0;
2919 :
2920 0 : Vector<Double> projparams(2);
2921 0 : projparams = 0.0;
2922 0 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2923 0 : Projection projTo( projection.type(), projparams );
2924 :
2925 : DirectionCoordinate
2926 0 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2927 : // projection,
2928 : projTo,
2929 0 : refCoord(0), refCoord(1),
2930 0 : deltas(0), deltas(1),
2931 : xform,
2932 0 : 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 0 : obslocation=obsPosition;
2940 0 : ObsInfo myobsinfo;
2941 0 : myobsinfo.setTelescope(telescop);
2942 0 : myobsinfo.setPointingCenter(mvPhaseCenter);
2943 0 : myobsinfo.setObsDate(obsEpoch);
2944 0 : 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 0 : if(spwids.nelements()==0)
2960 : {
2961 0 : Int nspw=msc.spectralWindow().nrow();
2962 0 : spwids.resize(nspw);
2963 0 : indgen(spwids);
2964 : }
2965 0 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2966 0 : Vector<Double> dataChanFreq, dataChanWidth;
2967 0 : std::vector<std::vector<Int> > averageWhichChan;
2968 0 : std::vector<std::vector<Int> > averageWhichSPW;
2969 0 : std::vector<std::vector<Double> > averageChanFrac;
2970 :
2971 0 : if(spwids.nelements()==1)
2972 : {
2973 0 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2974 0 : 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 0 : Double minDataFreq = min(dataChanFreq);
2985 0 : 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 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3037 0 : String cubemode;
3038 0 : if ( qrestfreq.getValue("Hz")==0 )
3039 : {
3040 0 : MSDopplerUtil msdoppler(msobj);
3041 0 : Vector<Double> restfreqvec;
3042 0 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3043 0 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3044 0 : 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 0 : }
3062 : Double refPix;
3063 0 : Vector<Double> chanFreq;
3064 0 : Vector<Double> chanFreqStep;
3065 0 : String specmode;
3066 :
3067 0 : 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 0 : 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 0 : Bool nonLinearFreq(false);
3083 0 : String veltype_p=veltype;
3084 0 : veltype_p.upcase();
3085 0 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3086 0 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3087 : {
3088 0 : nonLinearFreq= true;
3089 : }
3090 :
3091 0 : SpectralCoordinate mySpectral;
3092 : Double stepf;
3093 0 : if(!nonLinearFreq)
3094 : //TODO: velocity mode default start case (use last channels?)
3095 : {
3096 0 : Double startf=chanFreq[0];
3097 : //Double stepf=chanFreqStep[0];
3098 0 : if(chanFreq.nelements()==1)
3099 : {
3100 0 : stepf=chanFreqStep[0];
3101 : }
3102 : else
3103 : {
3104 0 : stepf=chanFreq[1]-chanFreq[0];
3105 : }
3106 0 : Double restf=qrestfreq.getValue("Hz");
3107 : //stepf=9e8;
3108 0 : 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 0 : 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 0 : 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 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3129 0 : 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 0 : uInt nrestfreq = restFreq.nelements();
3159 0 : 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 0 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3177 0 : if(whichStokes.nelements()==0)
3178 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3179 0 : StokesCoordinate myStokes(whichStokes);
3180 :
3181 : ////////////////// Build Full coordinate system.
3182 :
3183 : //CoordinateSystem csys;
3184 0 : csys.addCoordinate(myRaDec);
3185 0 : csys.addCoordinate(myStokes);
3186 0 : csys.addCoordinate(mySpectral);
3187 0 : csys.setObsInfo(myobsinfo);
3188 :
3189 : //store back csys to impars record
3190 : //cerr<<"save csys to csysRecord..."<<endl;
3191 0 : if(csysRecord.isDefined("coordsys"))
3192 0 : csysRecord.removeField("coordsys");
3193 0 : csys.save(csysRecord,"coordsys");
3194 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3195 : // imshape
3196 0 : imshape.resize(4);
3197 0 : imshape[0] = imsize[0];
3198 0 : imshape[1] = imsize[0];
3199 0 : imshape[2] = whichStokes.nelements();
3200 0 : 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 0 : } // end of else when coordsys record is not defined...
3234 :
3235 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3236 :
3237 0 : return csys;
3238 0 : }
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 0 : 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 0 : String inStart, inStep;
3531 0 : specmode = findSpecMode(mode);
3532 0 : String freqframe;
3533 0 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3534 0 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3535 :
3536 0 : refPix=0.0;
3537 0 : Bool descendingfreq(false);
3538 0 : Bool descendingoutfreq(false);
3539 :
3540 0 : 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 0 : else if ( mode=="mfs" ) {
3689 0 : chanFreq.resize(1);
3690 0 : chanFreqStep.resize(1);
3691 : //chanFreqStep[0] = freqmax - freqmin;
3692 0 : Double freqmean = (freqmin + freqmax)/2;
3693 0 : if (refFreq.getValue("Hz")==0) {
3694 0 : chanFreq[0] = freqmean;
3695 0 : refPix = 0.0;
3696 0 : 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 0 : if( nchan==-1 ) nchan=1;
3707 0 : 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 0 : return true;
3719 :
3720 0 : }//getImFreq
3721 : /////////////////////////
3722 0 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3723 0 : Double inStartFreq=-1.0;
3724 0 : String checkspecmode("");
3725 0 : if(mode.contains("cube")) {
3726 0 : checkspecmode = findSpecMode(mode);
3727 : }
3728 0 : 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 0 : return inStartFreq;
3748 :
3749 0 : }
3750 :
3751 0 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3752 : {
3753 0 : String specmode;
3754 0 : specmode="channel";
3755 0 : 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 0 : return specmode;
3770 0 : }
3771 :
3772 :
3773 0 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3774 : {
3775 0 : Vector<Int> whichStokes(0);
3776 0 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3777 0 : stokes=="RR" ||stokes=="LL" ||
3778 0 : stokes=="XX" || stokes=="YY" ) {
3779 0 : whichStokes.resize(1);
3780 0 : 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 0 : return whichStokes;
3810 0 : }// decidenpolplanes
3811 :
3812 0 : IPosition SynthesisParamsImage::shp() const
3813 : {
3814 0 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3815 :
3816 0 : 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 0 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3822 : }
3823 :
3824 0 : Record SynthesisParamsImage::getcsys() const
3825 : {
3826 0 : 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 0 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3883 : {
3884 0 : setDefaults();
3885 0 : }
3886 :
3887 0 : SynthesisParamsGrid::~SynthesisParamsGrid()
3888 : {
3889 0 : }
3890 :
3891 :
3892 0 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3893 : {
3894 0 : setDefaults();
3895 :
3896 0 : String err("");
3897 :
3898 : try {
3899 0 : err += readVal( inrec, String("imagename"), imageName );
3900 :
3901 : // FTMachine parameters
3902 0 : err += readVal( inrec, String("gridder"), gridder );
3903 0 : err += readVal( inrec, String("padding"), padding );
3904 0 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3905 0 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3906 0 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3907 0 : err += readVal( inrec, String("convfunc"), convFunc );
3908 :
3909 0 : err += readVal( inrec, String("vptable"), vpTable );
3910 :
3911 :
3912 :
3913 : // convert 'gridder' to 'ftmachine' and 'mtype'
3914 0 : ftmachine = "gridft";
3915 0 : mType = "default";
3916 0 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
3917 0 : ftmachine = "gridft";
3918 : }
3919 0 : else if ( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) &&
3920 0 : (wprojplanes>1 || wprojplanes==-1) ) {
3921 0 : 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 0 : String deconvolver;
3953 0 : err += readVal( inrec, String("deconvolver"), deconvolver );
3954 0 : if (deconvolver=="mtmfs") {
3955 0 : mType = "multiterm"; // Takes precedence over imagemosaic
3956 : }
3957 :
3958 : // facets
3959 0 : err += readVal( inrec, String("facets"), facets );
3960 : // chanchunks
3961 0 : err += readVal( inrec, String("chanchunks"), chanchunks );
3962 :
3963 : // Spectral interpolation
3964 0 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
3965 : // Track moving source ?
3966 0 : err += readVal( inrec, String("distance"), distance );
3967 0 : err += readVal( inrec, String("tracksource"), trackSource );
3968 0 : err += readVal( inrec, String("trackdir"), trackDir );
3969 :
3970 : // The extra params for WB-AWP
3971 0 : err += readVal( inrec, String("aterm"), aTermOn );
3972 0 : err += readVal( inrec, String("psterm"), psTermOn );
3973 0 : err += readVal( inrec, String("mterm"), mTermOn );
3974 0 : err += readVal( inrec, String("wbawp"), wbAWP );
3975 0 : err += readVal( inrec, String("cfcache"), cfCache );
3976 0 : err += readVal( inrec, String("usepointing"), usePointing );
3977 0 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3978 0 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3979 0 : err += readVal( inrec, String("conjbeams"), conjBeams );
3980 0 : err += readVal( inrec, String("computepastep"), computePAStep );
3981 0 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3982 :
3983 : // The extra params for single-dish
3984 0 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3985 0 : err += readVal( inrec, String("convertfirst"), convertFirst );
3986 0 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3987 0 : err += readVal( inrec, String("convsupport"), convSupport );
3988 0 : err += readVal( inrec, String("truncate"), truncateSize );
3989 0 : err += readVal( inrec, String("gwidth"), gwidth );
3990 0 : err += readVal( inrec, String("jwidth"), jwidth );
3991 0 : err += readVal( inrec, String("minweight"), minWeight );
3992 0 : 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 0 : if (ftmachine=="awprojectft" && cfCache=="") {
3998 0 : cfCache = imageName + ".cf";
3999 : }
4000 :
4001 0 : if ( ftmachine=="awprojectft" &&
4002 0 : 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 0 : err += verify();
4011 :
4012 0 : } catch(AipsError &x) {
4013 0 : err = err + x.getMesg() + "\n";
4014 0 : }
4015 :
4016 0 : if (err.length()>0) {
4017 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4018 : }
4019 :
4020 0 : }
4021 :
4022 0 : String SynthesisParamsGrid::verify() const
4023 : {
4024 0 : 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 0 : if ( imageName == "" ) {
4032 0 : err += "Please supply an image name\n";
4033 : }
4034 :
4035 0 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4036 0 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4037 0 : (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 0 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4049 0 : ( 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : return err;
4124 0 : }
4125 :
4126 0 : void SynthesisParamsGrid::setDefaults()
4127 : {
4128 0 : imageName="";
4129 : // FTMachine parameters
4130 : //ftmachine="GridFT";
4131 0 : ftmachine="gridft";
4132 0 : gridder=ftmachine;
4133 0 : padding=1.2;
4134 0 : useAutoCorr=false;
4135 0 : useDoublePrec=true;
4136 0 : wprojplanes=1;
4137 0 : convFunc="SF";
4138 0 : vpTable="";
4139 :
4140 : // facets
4141 0 : facets=1;
4142 :
4143 : // chanchunks
4144 0 : chanchunks=1;
4145 :
4146 : // Spectral Axis interpolation
4147 0 : interpolation=String("nearest");
4148 :
4149 : //mosaic use pointing
4150 0 : usePointing=false;
4151 : // Moving phase center ?
4152 0 : distance=Quantity(0,"m");
4153 0 : trackSource=false;
4154 0 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4155 :
4156 : // The extra params for WB-AWP
4157 0 : aTermOn = true;
4158 0 : psTermOn = true;
4159 0 : mTermOn = false;
4160 0 : wbAWP = true;
4161 0 : cfCache = "";
4162 0 : usePointing = false;
4163 0 : pointingOffsetSigDev.resize(0);
4164 : // pointingOffsetSigDev.set(30.0);
4165 0 : doPBCorr = true;
4166 0 : conjBeams = true;
4167 0 : computePAStep=360.0;
4168 0 : rotatePAStep=5.0;
4169 :
4170 : // extra params for single-dish
4171 0 : pointingDirCol = "";
4172 0 : convertFirst = "never";
4173 0 : skyPosThreshold = 0.0;
4174 0 : convSupport = -1;
4175 0 : truncateSize = Quantity(-1.0);
4176 0 : gwidth = Quantity(-1.0);
4177 0 : jwidth = Quantity(-1.0);
4178 0 : minWeight = 0.0;
4179 0 : clipMinMax = False;
4180 :
4181 : // Mapper type
4182 0 : mType = String("default");
4183 :
4184 0 : }
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 0 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4248 : {
4249 0 : setDefaults();
4250 0 : }
4251 :
4252 0 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4253 : {
4254 0 : }
4255 :
4256 :
4257 0 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4258 : {
4259 0 : setDefaults();
4260 :
4261 0 : String err("");
4262 :
4263 : try
4264 : {
4265 :
4266 0 : err += readVal( inrec, String("imagename"), imageName );
4267 0 : 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 0 : if( inrec.isDefined("startmodel") )
4273 : {
4274 0 : if( inrec.dataType("startmodel")==TpString )
4275 : {
4276 0 : String onemodel;
4277 0 : err += readVal( inrec, String("startmodel"), onemodel );
4278 0 : if( onemodel.length()>0 )
4279 : {
4280 0 : startModel.resize(1);
4281 0 : startModel[0] = onemodel;
4282 : }
4283 0 : else {startModel.resize();}
4284 0 : }
4285 0 : else if( inrec.dataType("startmodel")==TpArrayString ||
4286 0 : inrec.dataType("startmodel")==TpArrayBool)
4287 : {
4288 0 : 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 0 : err += readVal( inrec, String("id"), deconvolverId );
4297 0 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4298 :
4299 0 : err += readVal( inrec, String("scales"), scales );
4300 0 : err += readVal( inrec, String("scalebias"), scalebias );
4301 :
4302 0 : err += readVal( inrec, String("usemask"), maskType );
4303 0 : if( maskType=="auto-thresh" )
4304 : {
4305 0 : autoMaskAlgorithm = "thresh";
4306 : }
4307 0 : else if( maskType=="auto-thesh2" )
4308 : {
4309 0 : autoMaskAlgorithm = "thresh2";
4310 : }
4311 0 : else if( maskType=="auto-onebox" )
4312 : {
4313 0 : autoMaskAlgorithm = "onebox";
4314 : }
4315 0 : else if( maskType=="user" || maskType=="pb" )
4316 : {
4317 0 : autoMaskAlgorithm = "";
4318 : }
4319 :
4320 :
4321 0 : if( inrec.isDefined("mask") )
4322 : {
4323 0 : if( inrec.dataType("mask")==TpString )
4324 : {
4325 0 : 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 0 : if( inrec.isDefined("pbmask") )
4334 : {
4335 0 : err += readVal( inrec, String("pbmask"), pbMask );
4336 : }
4337 0 : if( inrec.isDefined("maskthreshold") )
4338 : {
4339 0 : if( inrec.dataType("maskthreshold")==TpString )
4340 : {
4341 0 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4342 : //deal with the case a string is a float value without unit
4343 0 : Quantity testThresholdString;
4344 0 : Quantity::read(testThresholdString,maskThreshold);
4345 0 : if( testThresholdString.getUnit()=="" )
4346 : {
4347 0 : if(testThresholdString.getValue()<1.0)
4348 : {
4349 0 : fracOfPeak = testThresholdString.getValue();
4350 0 : maskThreshold=String("");
4351 : }
4352 : }
4353 0 : }
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 0 : if( inrec.isDefined("maskresolution") )
4372 : {
4373 0 : if( inrec.dataType("maskresolution")==TpString )
4374 : {
4375 0 : err += readVal(inrec, String("maskresolution"), maskResolution );
4376 : //deal with the case a string is a float value without unit
4377 0 : Quantity testResolutionString;
4378 0 : Quantity::read(testResolutionString,maskResolution);
4379 0 : if( testResolutionString.getUnit()=="" )
4380 : {
4381 0 : maskResByBeam = testResolutionString.getValue();
4382 0 : maskResolution=String("");
4383 : }
4384 0 : }
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 0 : if( inrec.isDefined("nmask") )
4398 : {
4399 0 : if( inrec.dataType("nmask")==TpInt )
4400 : {
4401 0 : err+= readVal(inrec, String("nmask"), nMask );
4402 : }
4403 : else
4404 : {
4405 0 : err+= "nmask must be an integer\n";
4406 : }
4407 : }
4408 0 : if( inrec.isDefined("autoadjust") )
4409 : {
4410 0 : if( inrec.dataType("autoadjust")==TpBool )
4411 : {
4412 0 : 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 0 : if (inrec.isDefined("fusedthreshold"))
4421 : {
4422 0 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4423 0 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4424 : else
4425 0 : err += "fusedthreshold must be a float or double";
4426 : }
4427 0 : if (inrec.isDefined("specmode"))
4428 : {
4429 0 : if(inrec.dataType("specmode") == TpString)
4430 0 : 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 0 : if (inrec.isDefined("largestscale"))
4436 : {
4437 0 : if (inrec.dataType("largestscale") == TpInt)
4438 0 : 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 0 : if( inrec.isDefined("sidelobethreshold"))
4444 : {
4445 0 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4446 : {
4447 0 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4448 : }
4449 : else
4450 : {
4451 0 : err+= "sidelobethreshold must be a float or double";
4452 : }
4453 : }
4454 :
4455 0 : if( inrec.isDefined("noisethreshold"))
4456 : {
4457 0 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4458 : {
4459 0 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4460 : }
4461 : else
4462 : {
4463 0 : err+= "noisethreshold must be a float or double";
4464 : }
4465 : }
4466 0 : if( inrec.isDefined("lownoisethreshold"))
4467 : {
4468 0 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4469 : {
4470 0 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4471 : }
4472 : else
4473 : {
4474 0 : err+= "lownoisethreshold must be a float or double";
4475 : }
4476 : }
4477 0 : if( inrec.isDefined("negativethreshold"))
4478 : {
4479 0 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4480 : {
4481 0 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4482 : }
4483 : else
4484 : {
4485 0 : err+= "negativethreshold must be a float or double";
4486 : }
4487 : }
4488 0 : if( inrec.isDefined("smoothfactor"))
4489 : {
4490 0 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4491 : {
4492 0 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4493 : }
4494 : else
4495 : {
4496 0 : err+= "smoothfactor must be a float or double";
4497 : }
4498 : }
4499 0 : if( inrec.isDefined("minbeamfrac"))
4500 : {
4501 0 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4502 : {
4503 0 : 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 0 : if( inrec.isDefined("cutthreshold"))
4517 : {
4518 0 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4519 : {
4520 0 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4521 : }
4522 : else {
4523 0 : err+= "cutthreshold must be a float or double";
4524 : }
4525 : }
4526 0 : if( inrec.isDefined("growiterations"))
4527 : {
4528 0 : if (inrec.dataType("growiterations")==TpInt) {
4529 0 : err+= readVal(inrec, String("growiterations"), growIterations );
4530 : }
4531 : else {
4532 0 : err+= "growiterations must be an integer\n";
4533 : }
4534 : }
4535 0 : if( inrec.isDefined("dogrowprune"))
4536 : {
4537 0 : if (inrec.dataType("dogrowprune")==TpBool) {
4538 0 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4539 : }
4540 : else {
4541 0 : err+= "dogrowprune must be a bool\n";
4542 : }
4543 : }
4544 0 : if( inrec.isDefined("minpercentchange"))
4545 : {
4546 0 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4547 0 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4548 : }
4549 : else {
4550 0 : err+= "minpercentchange must be a float or double";
4551 : }
4552 : }
4553 0 : if( inrec.isDefined("verbose"))
4554 : {
4555 0 : if (inrec.dataType("verbose")==TpBool ) {
4556 0 : err+= readVal(inrec, String("verbose"), verbose);
4557 : }
4558 : else {
4559 0 : err+= "verbose must be a bool";
4560 : }
4561 : }
4562 0 : if( inrec.isDefined("fastnoise"))
4563 : {
4564 0 : if (inrec.dataType("fastnoise")==TpBool ) {
4565 0 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4566 : }
4567 : else {
4568 0 : err+= "fastnoise must be a bool";
4569 : }
4570 : }
4571 0 : if( inrec.isDefined("nsigma") )
4572 : {
4573 0 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4574 0 : 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 0 : if( inrec.isDefined("noRequireSumwt") )
4587 : {
4588 0 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4589 0 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4590 : }
4591 : else {
4592 0 : err+= "noRequireSumwt must be a bool";
4593 : }
4594 : }
4595 0 : if( inrec.isDefined("fullsummary") )
4596 : {
4597 0 : if (inrec.dataType("fullsummary")==TpBool) {
4598 0 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4599 : }
4600 : else {
4601 0 : err+= "fullsummary must be a bool";
4602 : }
4603 : }
4604 0 : if( inrec.isDefined("restoringbeam") )
4605 : {
4606 0 : String errinfo("");
4607 : try {
4608 :
4609 0 : 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 0 : 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 0 : }// if isdefined(restoringbeam)
4658 :
4659 0 : if( inrec.isDefined("interactive") )
4660 : {
4661 0 : if( inrec.dataType("interactive")==TpBool )
4662 0 : {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 0 : err += verify();
4670 :
4671 : }
4672 0 : catch(AipsError &x)
4673 : {
4674 0 : err = err + x.getMesg() + "\n";
4675 0 : }
4676 :
4677 0 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4678 :
4679 0 : }
4680 :
4681 0 : String SynthesisParamsDeconv::verify() const
4682 : {
4683 0 : String err;
4684 :
4685 0 : 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 0 : 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 0 : if( pbMask >= 1.0)
4695 0 : {err += "pbmask must be < 1.0 \n"; }
4696 0 : else if( pbMask < 0.0)
4697 0 : {err += "pbmask must be a positive value \n"; }
4698 :
4699 0 : 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 0 : if ( fracOfPeak >= 1.0)
4708 0 : {err += "fracofpeak must be < 1.0 \n"; }
4709 0 : else if ( fracOfPeak < 0.0)
4710 0 : {err += "fracofpeak must be a positive value \n"; }
4711 :
4712 0 : return err;
4713 0 : }
4714 :
4715 0 : void SynthesisParamsDeconv::setDefaults()
4716 : {
4717 0 : imageName="";
4718 0 : algorithm="hogbom";
4719 0 : startModel=Vector<String>(0);
4720 0 : deconvolverId=0;
4721 0 : nTaylorTerms=1;
4722 0 : scales.resize(1); scales[0]=0.0;
4723 0 : scalebias=0.6;
4724 0 : maskType="none";
4725 0 : maskString="";
4726 0 : maskList.resize(1); maskList[0]="";
4727 0 : pbMask=0.0;
4728 0 : autoMaskAlgorithm="thresh";
4729 0 : maskThreshold="";
4730 0 : maskResolution="";
4731 0 : fracOfPeak=0.0;
4732 0 : nMask=0;
4733 0 : interactive=false;
4734 0 : autoAdjust=False;
4735 0 : fusedThreshold = 0.0;
4736 0 : specmode="mfs";
4737 0 : largestscale = -1;
4738 0 : }
4739 :
4740 0 : Record SynthesisParamsDeconv::toRecord() const
4741 : {
4742 0 : Record decpar;
4743 :
4744 0 : decpar.define("imagename", imageName);
4745 0 : decpar.define("deconvolver", algorithm);
4746 0 : decpar.define("startmodel",startModel);
4747 0 : decpar.define("id",deconvolverId);
4748 0 : decpar.define("nterms",nTaylorTerms);
4749 0 : decpar.define("scales",scales);
4750 0 : decpar.define("scalebias",scalebias);
4751 0 : decpar.define("usemask",maskType);
4752 0 : decpar.define("fusedthreshold", fusedThreshold);
4753 0 : decpar.define("specmode", specmode);
4754 0 : decpar.define("largestscale", largestscale);
4755 0 : if( maskList.nelements()==1 && maskList[0]=="")
4756 : {
4757 0 : decpar.define("mask",maskString);
4758 : }
4759 : else {
4760 0 : decpar.define("mask",maskList);
4761 : }
4762 0 : decpar.define("pbmask",pbMask);
4763 0 : if (fracOfPeak > 0.0)
4764 : {
4765 0 : decpar.define("maskthreshold",fracOfPeak);
4766 : }
4767 : else
4768 : {
4769 0 : decpar.define("maskthreshold",maskThreshold);
4770 : }
4771 0 : decpar.define("maskresolution",maskResolution);
4772 0 : decpar.define("nmask",nMask);
4773 0 : decpar.define("autoadjust",autoAdjust);
4774 0 : decpar.define("sidelobethreshold",sidelobeThreshold);
4775 0 : decpar.define("noisethreshold",noiseThreshold);
4776 0 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4777 0 : decpar.define("negativethreshold",negativeThreshold);
4778 0 : decpar.define("smoothfactor",smoothFactor);
4779 0 : decpar.define("minbeamfrac",minBeamFrac);
4780 0 : decpar.define("cutthreshold",cutThreshold);
4781 0 : decpar.define("growiterations",growIterations);
4782 0 : decpar.define("dogrowprune",doGrowPrune);
4783 0 : decpar.define("minpercentchange",minPercentChange);
4784 0 : decpar.define("verbose", verbose);
4785 0 : decpar.define("fastnoise", fastnoise);
4786 0 : decpar.define("interactive",interactive);
4787 0 : decpar.define("nsigma",nsigma);
4788 0 : decpar.define("noRequireSumwt",noRequireSumwt);
4789 0 : decpar.define("fullsummary",fullsummary);
4790 :
4791 0 : return decpar;
4792 0 : }
4793 :
4794 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4795 :
4796 :
4797 : } //# NAMESPACE CASA - END
4798 :
|
