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 20 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 20 : }
89 :
90 20 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 20 : }
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 5958 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 5958 : Int N=vb.nRows(),M=-1;
108 47454 : for(Int i=0;i<N;i++)
109 : {
110 47454 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 5958 : {M++;break;}
112 : }
113 5958 : 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 70 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 70 : Int n=npix;
121 :
122 70 : if( n%2 !=0 ){ n+= 1; }
123 :
124 70 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 380 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 310 : if( fac[k]>5 )
129 : {
130 10 : val = fac[k];
131 20 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 10 : fac[k] = val;
133 : }
134 : }
135 70 : newlarge=product(fac);
136 80 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 10 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 70 : return newlarge;
141 70 : }
142 :
143 : // Return the prime factors of the given number
144 100 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 100 : Vector<uInt> factors;
147 :
148 100 : Int lastresult = n;
149 100 : Int sqlast = int(sqrt(n))+1;
150 :
151 100 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 100 : Int c=2;
153 : while(1)
154 : {
155 470 : if( lastresult == 1 || c > sqlast ) { break; }
156 370 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 610 : if( c>sqlast){ c=lastresult; break; }
160 512 : if( lastresult % c == 0 ) { break; }
161 240 : c += 1;
162 : }
163 370 : factors.resize( factors.nelements()+1, true );
164 370 : factors[ factors.nelements()-1 ] = c;
165 370 : lastresult /= c;
166 : }
167 100 : 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 100 : return factors;
185 0 : }
186 :
187 :
188 3 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 6 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 3 : 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 3 : std::shared_ptr<SIImageStore> imstore;
199 3 : if( nterms>1 )
200 1 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 2 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 3 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 3 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 3 : imstore->printBeamSet();
210 :
211 3 : imstore->releaseLocks();
212 :
213 3 : return true;
214 3 : }
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 28 : 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 28 : Record outRec;
705 28 : outRec.define("type", type);
706 28 : outRec.define("rmode", rmode);
707 28 : Record quantRec;
708 28 : QuantumHolder(noise).toRecord(quantRec);
709 28 : outRec.defineRecord("noise", quantRec);
710 28 : outRec.define("robust", robust);
711 28 : QuantumHolder(fieldofview).toRecord(quantRec);
712 28 : outRec.defineRecord("fieldofview", quantRec);
713 28 : outRec.define("npixels", npixels);
714 28 : outRec.define("multifield", multiField);
715 28 : outRec.define("usecubebriggs", useCubeBriggs);
716 28 : outRec.define("filtertype", filtertype);
717 28 : QuantumHolder(filterbmaj).toRecord(quantRec);
718 28 : outRec.defineRecord("filterbmaj", quantRec);
719 28 : QuantumHolder(filterbmin).toRecord(quantRec);
720 28 : outRec.defineRecord("filterbmin", quantRec);
721 28 : QuantumHolder(filterbpa).toRecord(quantRec);
722 28 : outRec.defineRecord("filterbpa", quantRec);
723 28 : outRec.define("fracBW", fracBW);
724 :
725 56 : return outRec;
726 28 : }
727 21 : 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 21 : QuantumHolder qh;
734 21 : String err;
735 21 : if(!inRec.isDefined("type"))
736 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
737 21 : inRec.get("type", type);
738 21 : inRec.get("rmode", rmode);
739 21 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
740 0 : throw(AipsError("Error in reading noise param"));
741 21 : noise=qh.asQuantity();
742 21 : inRec.get("robust", robust);
743 21 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
744 0 : throw(AipsError("Error in reading fieldofview param"));
745 21 : fieldofview=qh.asQuantity();
746 21 : inRec.get("npixels", npixels);
747 21 : inRec.get("multifield", multiField);
748 21 : inRec.get("usecubebriggs", useCubeBriggs);
749 21 : inRec.get("filtertype", filtertype);
750 21 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
751 0 : throw(AipsError("Error in reading filterbmaj param"));
752 21 : filterbmaj=qh.asQuantity();
753 21 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
754 0 : throw(AipsError("Error in reading filterbmin param"));
755 21 : filterbmin=qh.asQuantity();
756 21 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
757 0 : throw(AipsError("Error in reading filterbpa param"));
758 21 : filterbpa=qh.asQuantity();
759 21 : inRec.get("fracBW", fracBW);
760 :
761 :
762 :
763 21 : }
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 338 : 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 338 : Bool isOn = false;
924 338 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
925 338 : if (!isOn)
926 338 : 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 7 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
1398 : {
1399 7 : MVAngle mvRa=direction.getAngle().getValue()(0);
1400 7 : MVAngle mvDec=direction.getAngle().getValue()(1);
1401 7 : ostringstream oos;
1402 7 : oos << " ";
1403 7 : Int widthRA=20;
1404 7 : Int widthDec=20;
1405 7 : oos.setf(ios::left, ios::adjustfield);
1406 7 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
1407 7 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
1408 7 : oos << " "
1409 7 : << MDirection::showType(direction.getRefPtr()->getType());
1410 14 : return String(oos);
1411 7 : }
1412 :
1413 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1414 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1415 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
1416 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1417 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1418 :
1419 : // Read String from Record
1420 1608 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1421 : {
1422 1608 : if( rec.isDefined( id ) )
1423 : {
1424 1496 : String inval("");
1425 1496 : if( rec.dataType( id )==TpString )
1426 1496 : { 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 1496 : if(inval.length()>0){val=inval;}
1432 1496 : return String("");
1433 : }
1434 0 : else { return String(id + " must be a string\n"); }
1435 1496 : }
1436 112 : else { return String("");}
1437 : }
1438 :
1439 : // Read Integer from Record
1440 411 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1441 : {
1442 411 : if( rec.isDefined( id ) )
1443 : {
1444 400 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1445 0 : else { return String(id + " must be an integer\n"); }
1446 : }
1447 11 : else { return String(""); }
1448 : }
1449 :
1450 : // Read Float from Record
1451 547 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1452 : {
1453 547 : if( rec.isDefined( id ) )
1454 : {
1455 498 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1456 498 : { rec.get( id , val ); return String(""); }
1457 0 : else { return String(id + " must be a float\n"); }
1458 : }
1459 49 : else { return String(""); }
1460 : }
1461 :
1462 : // Read Bool from Record
1463 809 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1464 : {
1465 809 : if( rec.isDefined( id ) )
1466 : {
1467 711 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1468 0 : else { return String(id + " must be a bool\n"); }
1469 : }
1470 98 : else{ return String(""); }
1471 : }
1472 :
1473 : // Read Vector<Int> from Record
1474 21 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
1475 : {
1476 21 : if( rec.isDefined( id ) )
1477 : {
1478 21 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1479 21 : { 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 66 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1493 : {
1494 66 : if( rec.isDefined( id ) )
1495 : {
1496 66 : if( rec.dataType( id )==TpArrayFloat )
1497 : {
1498 47 : 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 19 : else if ( rec.dataType( id ) ==TpArrayDouble )
1512 : {
1513 7 : Vector<Double> vec; rec.get( id, vec );
1514 7 : val.resize(vec.nelements());
1515 21 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1516 7 : return String("");
1517 7 : }
1518 12 : else if ( rec.dataType( id ) ==TpArrayInt )
1519 : {
1520 5 : Vector<Int> vec; rec.get( id, vec );
1521 5 : val.resize(vec.nelements());
1522 28 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1523 5 : return String("");
1524 5 : }
1525 7 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1526 : {
1527 7 : Vector<Bool> vec; rec.get( id, vec );
1528 7 : 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 7 : }
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 82 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1539 : {
1540 82 : if( rec.isDefined( id ) )
1541 : {
1542 82 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1543 82 : { 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 219 : 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 219 : 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 53 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1570 : {
1571 : try
1572 : {
1573 53 : String tmpRF, tmpRA, tmpDEC;
1574 53 : std::istringstream iss(instr);
1575 53 : iss >> tmpRF >> tmpRA >> tmpDEC;
1576 53 : 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 53 : casacore::Quantity tmpQRA;
1583 53 : casacore::Quantity tmpQDEC;
1584 53 : casacore::Quantity::read(tmpQRA, tmpRA);
1585 53 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1586 :
1587 53 : 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 53 : Bool status = MDirection::getType(theRF, tmpRF);
1598 53 : if (!status) {
1599 0 : throw AipsError();
1600 : }
1601 53 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1602 53 : return String("");
1603 53 : }
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 175 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1612 : {
1613 175 : if( rec.isDefined( id ) )
1614 : {
1615 161 : if( rec.dataType( id )==TpString )
1616 161 : { 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 14 : else{ return String(""); }
1620 : }
1621 :
1622 : // Read MDirection from a Record string
1623 67 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1624 : {
1625 67 : if( rec.isDefined( id ) )
1626 : {
1627 53 : if( rec.dataType( id )==TpString )
1628 53 : { 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 14 : else{ return String(""); }
1632 : }
1633 :
1634 : // Convert MDirection to String
1635 42 : String SynthesisParams::MDirectionToString(MDirection val) const
1636 : {
1637 42 : MVDirection mvpc( val.getAngle() );
1638 42 : MVAngle ra = mvpc.get()(0);
1639 42 : MVAngle dec = mvpc.get()(1);
1640 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1641 126 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1642 42 : }
1643 :
1644 : // Convert Quantity to String
1645 139 : String SynthesisParams::QuantityToString(Quantity val) const
1646 : {
1647 139 : 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 139 : ss.precision(std::numeric_limits<double>::digits10+2);
1652 139 : ss << val;
1653 417 : 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 139 : }
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 84 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1683 : {
1684 84 : setDefaults();
1685 84 : }
1686 :
1687 84 : SynthesisParamsSelect::~SynthesisParamsSelect()
1688 : {
1689 84 : }
1690 :
1691 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1692 0 : operator=(other);
1693 0 : }
1694 :
1695 35 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1696 35 : if(this!=&other) {
1697 35 : msname=other.msname;
1698 35 : spw=other.spw;
1699 35 : freqbeg=other.freqbeg;
1700 35 : freqend=other.freqend;
1701 35 : freqframe=other.freqframe;
1702 35 : field=other.field;
1703 35 : antenna=other.antenna;
1704 35 : timestr=other.timestr;
1705 35 : scan=other.scan;
1706 35 : obs=other.obs;
1707 35 : state=other.state;
1708 35 : uvdist=other.uvdist;
1709 35 : taql=other.taql;
1710 35 : usescratch=other.usescratch;
1711 35 : readonly=other.readonly;
1712 35 : incrmodel=other.incrmodel;
1713 35 : datacolumn=other.datacolumn;
1714 :
1715 : }
1716 35 : return *this;
1717 : }
1718 :
1719 49 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1720 : {
1721 49 : setDefaults();
1722 :
1723 49 : String err("");
1724 :
1725 : try
1726 : {
1727 :
1728 49 : err += readVal( inrec, String("msname"), msname );
1729 :
1730 49 : err += readVal( inrec, String("readonly"), readonly );
1731 49 : err += readVal( inrec, String("usescratch"), usescratch );
1732 :
1733 : // override with entries from savemodel.
1734 49 : String savemodel("");
1735 49 : err += readVal( inrec, String("savemodel"), savemodel );
1736 49 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1737 35 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1738 35 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1739 :
1740 49 : err += readVal( inrec, String("incrmodel"), incrmodel );
1741 :
1742 49 : err += readVal( inrec, String("spw"), spw );
1743 :
1744 : /// -------------------------------------------------------------------------------------
1745 : // Why are these params here ? Repeats of defineimage.
1746 49 : err += readVal( inrec, String("freqbeg"), freqbeg);
1747 49 : err += readVal( inrec, String("freqend"), freqend);
1748 :
1749 49 : String freqframestr( MFrequency::showType(freqframe) );
1750 49 : err += readVal( inrec, String("outframe"), freqframestr);
1751 49 : if( ! MFrequency::getType(freqframe, freqframestr) )
1752 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1753 : /// -------------------------------------------------------------------------------------
1754 :
1755 49 : err += readVal( inrec, String("field"),field);
1756 49 : err += readVal( inrec, String("antenna"),antenna);
1757 49 : err += readVal( inrec, String("timestr"),timestr);
1758 49 : err += readVal( inrec, String("scan"),scan);
1759 49 : err += readVal( inrec, String("obs"),obs);
1760 49 : err += readVal( inrec, String("state"),state);
1761 49 : err += readVal( inrec, String("uvdist"),uvdist);
1762 49 : err += readVal( inrec, String("taql"),taql);
1763 :
1764 49 : err += readVal( inrec, String("datacolumn"),datacolumn);
1765 :
1766 49 : err += verify();
1767 :
1768 49 : }
1769 0 : catch(AipsError &x)
1770 : {
1771 0 : err = err + x.getMesg() + "\n";
1772 0 : }
1773 :
1774 49 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1775 :
1776 49 : }
1777 :
1778 49 : String SynthesisParamsSelect::verify() const
1779 : {
1780 49 : String err;
1781 : // Does the MS exist on disk.
1782 49 : Directory thems( msname );
1783 49 : if( thems.exists() )
1784 : {
1785 : // Is it readable ?
1786 49 : if( ! thems.isReadable() )
1787 0 : { err += "MS " + msname + " is not readable.\n"; }
1788 : // Depending on 'readonly', is the MS writable ?
1789 49 : 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 98 : return err;
1796 49 : }
1797 :
1798 133 : void SynthesisParamsSelect::setDefaults()
1799 : {
1800 133 : msname="";
1801 133 : spw="";
1802 133 : freqbeg="";
1803 133 : freqend="";
1804 133 : MFrequency::getType(freqframe,"LSRK");
1805 133 : field="";
1806 133 : antenna="";
1807 133 : timestr="";
1808 133 : scan="";
1809 133 : obs="";
1810 133 : state="";
1811 133 : uvdist="";
1812 133 : taql="";
1813 133 : usescratch=false;
1814 133 : readonly=true;
1815 133 : incrmodel=false;
1816 133 : datacolumn="corrected";
1817 133 : }
1818 :
1819 35 : Record SynthesisParamsSelect::toRecord()const
1820 : {
1821 35 : Record selpar;
1822 35 : selpar.define("msname",msname);
1823 35 : selpar.define("spw",spw);
1824 35 : selpar.define("freqbeg",freqbeg);
1825 35 : selpar.define("freqend",freqend);
1826 35 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1827 : //looks like fromRecord looks for outframe !
1828 35 : selpar.define("outframe", MFrequency::showType(freqframe));
1829 35 : selpar.define("field",field);
1830 35 : selpar.define("antenna",antenna);
1831 35 : selpar.define("timestr",timestr);
1832 35 : selpar.define("scan",scan);
1833 35 : selpar.define("obs",obs);
1834 35 : selpar.define("state",state);
1835 35 : selpar.define("uvdist",uvdist);
1836 35 : selpar.define("taql",taql);
1837 35 : selpar.define("usescratch",usescratch);
1838 35 : selpar.define("readonly",readonly);
1839 35 : selpar.define("incrmodel",incrmodel);
1840 35 : selpar.define("datacolumn",datacolumn);
1841 :
1842 35 : return selpar;
1843 0 : }
1844 :
1845 :
1846 : /////////////////////// Image Parameters
1847 :
1848 110 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1849 : {
1850 110 : setDefaults();
1851 110 : }
1852 :
1853 110 : SynthesisParamsImage::~SynthesisParamsImage()
1854 : {
1855 110 : }
1856 :
1857 70 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1858 70 : if(this != &other){
1859 70 : imageName=other.imageName;
1860 70 : stokes=other.stokes;
1861 70 : startModel.resize(); startModel=other.startModel;
1862 70 : imsize.resize(); imsize=other.imsize;
1863 70 : cellsize.resize(); cellsize=other.cellsize;
1864 70 : projection=other.projection;
1865 70 : useNCP=other.useNCP;
1866 70 : phaseCenter=other.phaseCenter;
1867 70 : phaseCenterFieldId=other.phaseCenterFieldId;
1868 70 : obslocation=other.obslocation;
1869 70 : pseudoi=other.pseudoi;
1870 70 : nchan=other.nchan;
1871 70 : nTaylorTerms=other.nTaylorTerms;
1872 70 : chanStart=other.chanStart;
1873 70 : chanStep=other.chanStep;
1874 70 : freqStart=other.freqStart;
1875 70 : freqStep=other.freqStep;
1876 70 : refFreq=other.refFreq;
1877 70 : velStart=other.velStart;
1878 70 : velStep=other.velStep;
1879 70 : freqFrame=other.freqFrame;
1880 70 : mFreqStart=other.mFreqStart;
1881 70 : mFreqStep=other.mFreqStep;
1882 70 : mVelStart=other.mVelStart;
1883 70 : mVelStep=other.mVelStep;
1884 70 : restFreq.resize(); restFreq=other.restFreq;
1885 70 : start=other.start;
1886 70 : step=other.step;
1887 70 : frame=other.frame;
1888 70 : veltype=other.veltype;
1889 70 : mode=other.mode;
1890 70 : reffreq=other.reffreq;
1891 70 : sysvel=other.sysvel;
1892 70 : sysvelframe=other.sysvelframe;
1893 70 : sysvelvalue=other.sysvelvalue;
1894 70 : qmframe=other.qmframe;
1895 70 : mveltype=other.mveltype;
1896 70 : tststr=other.tststr;
1897 70 : startRecord=other.startRecord;
1898 70 : stepRecord=other.stepRecord;
1899 70 : reffreqRecord=other.reffreqRecord;
1900 70 : sysvelRecord=other.sysvelRecord;
1901 70 : restfreqRecord=other.restfreqRecord;
1902 70 : csysRecord=other.csysRecord;
1903 70 : csys=other.csys;
1904 70 : imshape.resize(); imshape=other.imshape;
1905 70 : freqFrameValid=other.freqFrameValid;
1906 70 : overwrite=other.overwrite;
1907 70 : deconvolver=other.deconvolver;
1908 70 : distance=other.distance;
1909 70 : trackDir=other.trackDir;
1910 70 : trackSource=other.trackSource;
1911 70 : movingSource=other.movingSource;
1912 :
1913 :
1914 :
1915 : }
1916 :
1917 70 : return *this;
1918 :
1919 : }
1920 :
1921 28 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec) {
1922 28 : fromRecord(inrec, False);
1923 28 : }
1924 :
1925 35 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec,
1926 : const casacore::Bool isSingleDish) {
1927 35 : setDefaults();
1928 35 : String err("");
1929 70 : LogIO os(LogOrigin("SynthesisParamsImage", "fromRecord"));
1930 :
1931 : try {
1932 35 : err += readVal(inrec, String("imagename"), imageName);
1933 :
1934 : //// imsize
1935 35 : if (inrec.isDefined("imsize")) {
1936 35 : DataType tp = inrec.dataType("imsize");
1937 35 : if (tp == TpInt) { // A single integer for both dimensions.
1938 : Int npix;
1939 7 : inrec.get("imsize", npix);
1940 7 : imsize.resize(2);
1941 7 : imsize.set(npix);
1942 28 : } else if (tp == TpArrayInt) { // An integer array : [ nx ] or [ nx, ny ]
1943 28 : Vector<Int> ims;
1944 28 : inrec.get("imsize", ims);
1945 28 : if (ims.nelements() == 1) { // [ nx ]
1946 0 : imsize.set(ims[0]);
1947 28 : } else if (ims.nelements() == 2) { // [ nx, ny ]
1948 28 : imsize[0] = ims[0];
1949 28 : 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 28 : } 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 35 : if (inrec.isDefined("cell")) {
1960 : try {
1961 35 : DataType tp = inrec.dataType("cell");
1962 35 : if (tp == TpInt ||
1963 35 : tp == TpFloat ||
1964 : tp == TpDouble) {
1965 0 : Double cell = inrec.asDouble("cell");
1966 0 : cellsize.set(Quantity(cell, "arcsec"));
1967 35 : } else if (tp == TpArrayInt ||
1968 35 : 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 35 : } else if (tp == TpString) {
1981 12 : String cell;
1982 12 : inrec.get("cell", cell);
1983 12 : Quantity qcell;
1984 12 : err += stringToQuantity(cell, qcell);
1985 12 : cellsize.set(qcell);
1986 35 : } else if (tp == TpArrayString) {
1987 23 : Array<String> cells;
1988 23 : inrec.get("cell", cells);
1989 23 : Vector<String> vcells(cells);
1990 23 : if (cells.nelements() == 1) { // [ cellx ]
1991 0 : Quantity qcell;
1992 0 : err += stringToQuantity(vcells[0], qcell);
1993 0 : cellsize.set(qcell);
1994 23 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1995 23 : err += stringToQuantity(vcells[0], cellsize[0]);
1996 23 : 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 23 : } 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 35 : err += readVal(inrec, String("stokes"), stokes);
2010 35 : if (stokes.matches("pseudoI")) {
2011 1 : stokes = "I";
2012 1 : pseudoi = true;
2013 : } else {
2014 34 : pseudoi = false;
2015 : }
2016 :
2017 : /// PseudoI
2018 :
2019 : ////nchan
2020 35 : err += readVal(inrec, String("nchan"), nchan);
2021 :
2022 : /// phaseCenter (as a string) . // Add INT support later.
2023 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2024 35 : if (inrec.isDefined("phasecenter")) {
2025 35 : String pcent("");
2026 35 : if (inrec.dataType("phasecenter") == TpString) {
2027 35 : inrec.get("phasecenter", pcent);
2028 35 : if (pcent.length() > 0) { // if it's zero length, it means 'figure out from first field in MS'.
2029 32 : err += readVal(inrec, String("phasecenter"), phaseCenter);
2030 32 : 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 32 : err += readVal(inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2034 : } else {
2035 3 : phaseCenterFieldId = -2; // deal with this later in buildCoordinateSystem to assign the first selected field
2036 : }
2037 : }
2038 35 : if (inrec.dataType("phasecenter") == TpInt ||
2039 105 : inrec.dataType("phasecenter") == TpFloat ||
2040 70 : 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 35 : if ((inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter") != TpInt
2045 35 : && inrec.dataType("phasecenter") != TpFloat && inrec.dataType("phasecenter") != TpDouble)) {
2046 0 : err += String("Cannot set phasecenter. Please specify a string or int\n");
2047 : }
2048 35 : }
2049 :
2050 : //// Projection
2051 35 : if (inrec.isDefined("projection")) {
2052 35 : if (inrec.dataType("projection") == TpString) {
2053 35 : String pstr;
2054 35 : inrec.get("projection", pstr);
2055 : try {
2056 35 : if (pstr.matches("NCP")) {
2057 0 : pstr = "SIN";
2058 0 : useNCP = true;
2059 : }
2060 35 : projection = Projection::type(pstr);
2061 0 : } catch (AipsError &x) {
2062 0 : err += String("Invalid projection code : " + pstr);
2063 0 : }
2064 35 : } else {
2065 0 : err += "projection must be a string\n";
2066 : }
2067 : } //projection
2068 :
2069 : // Frequency frame stuff.
2070 35 : err += readVal(inrec, String("specmode"), mode);
2071 : // Alias for 'mfs' is 'cont'
2072 35 : if (mode == "cont") mode = "mfs";
2073 :
2074 35 : err += readVal(inrec, String("outframe"), frame);
2075 35 : qmframe = "";
2076 : // mveltype is only set when start/step is given in mdoppler
2077 35 : mveltype = "";
2078 : //start
2079 35 : String startType("");
2080 35 : String widthType("");
2081 35 : if (inrec.isDefined("start")) {
2082 35 : if (inrec.dataType("start") == TpInt) {
2083 7 : err += readVal(inrec, String("start"), chanStart);
2084 7 : start = String::toString(chanStart);
2085 7 : startType = "chan";
2086 28 : } else if (inrec.dataType("start") == TpString) {
2087 28 : err += readVal(inrec, String("start"), start);
2088 28 : if (start.contains("Hz")) {
2089 0 : stringToQuantity(start, freqStart);
2090 0 : startType = "freq";
2091 28 : } 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 35 : if (inrec.isDefined("width")) {
2152 35 : if (inrec.dataType("width") == TpInt) {
2153 7 : err += readVal(inrec, String("width"), chanStep);
2154 7 : step = String::toString(chanStep);
2155 7 : widthType = "chan";
2156 28 : } else if (inrec.dataType("width") == TpString) {
2157 28 : err += readVal(inrec, String("width"), step);
2158 28 : if (step.contains("Hz")) {
2159 0 : stringToQuantity(step, freqStep);
2160 0 : widthType = "freq";
2161 28 : } 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 35 : 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 35 : if (inrec.isDefined("reffreq")) {
2226 35 : if (inrec.dataType("reffreq") == TpString) {
2227 35 : 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 35 : err += readVal(inrec, String("veltype"), veltype);
2254 35 : veltype = mveltype != "" ? mveltype : veltype;
2255 : // sysvel (String, Quantity)
2256 35 : if (inrec.isDefined("sysvel")) {
2257 35 : if (inrec.dataType("sysvel") == TpString) {
2258 35 : 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 35 : err += readVal(inrec, String("sysvelframe"), sysvelframe);
2270 :
2271 : // rest frequencies (record or vector of Strings)
2272 35 : if (inrec.isDefined("restfreq")) {
2273 35 : Vector<String> rfreqs(0);
2274 35 : Record restfreqSubRecord;
2275 35 : 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 35 : } else if (inrec.dataType("restfreq") == TpArrayString) {
2286 21 : err += readVal(inrec, String("restfreq"), rfreqs);
2287 14 : } else if (inrec.dataType("restfreq") == TpString) {
2288 0 : rfreqs.resize(1);
2289 0 : err += readVal(inrec, String("restfreq"), rfreqs[0]);
2290 : }
2291 35 : restFreq.resize(rfreqs.nelements());
2292 35 : for (uInt fr = 0; fr < rfreqs.nelements(); fr++) {
2293 0 : err += stringToQuantity(rfreqs[fr], restFreq[fr]);
2294 : }
2295 35 : } // if def restfreq
2296 :
2297 : // optional - coordsys, imshape
2298 : // if exist use them. May need a consistency check with the rest of impars?
2299 35 : if (inrec.isDefined("csys")) {
2300 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
2301 21 : if (inrec.dataType("csys") == TpRecord) {
2302 : //csysRecord = inrec.subRecord("csys");
2303 21 : csysRecord.defineRecord("coordsys", inrec.subRecord("csys"));
2304 : }
2305 21 : if (inrec.isDefined("imshape")) {
2306 21 : if (inrec.dataType("imshape") == TpArrayInt) {
2307 21 : 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 35 : String freqframestr = (frame != "" && qmframe != "") ? qmframe : frame;
2318 35 : if (frame != "" && !MFrequency::getType(freqFrame, freqframestr)) {
2319 0 : err += "Invalid Frequency Frame " + freqframestr;
2320 : }
2321 35 : err += readVal(inrec, String("restart"), overwrite);
2322 :
2323 35 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2324 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2325 35 : if (inrec.isDefined("startmodel")) {
2326 35 : if (inrec.dataType("startmodel") == TpString) {
2327 14 : String onemodel;
2328 14 : err += readVal(inrec, String("startmodel"), onemodel);
2329 14 : if (onemodel.length() > 0) {
2330 0 : startModel.resize(1);
2331 0 : startModel[0] = onemodel;
2332 : } else {
2333 14 : startModel.resize();
2334 : }
2335 56 : } else if (inrec.dataType("startmodel") == TpArrayString ||
2336 21 : inrec.dataType("startmodel") == TpArrayBool) {
2337 21 : 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 35 : err += readVal(inrec, String("nterms"), nTaylorTerms);
2344 35 : err += readVal(inrec, String("deconvolver"), deconvolver);
2345 :
2346 : // Force nchan=1 for anything other than cube modes...
2347 35 : if (mode == "mfs") nchan = 1;
2348 : //read obslocation
2349 35 : if (inrec.isDefined("obslocation_rec")) {
2350 21 : String errorobs;
2351 21 : const Record obsrec = inrec.asRecord("obslocation_rec");
2352 21 : MeasureHolder mh;
2353 21 : if (!mh.fromRecord(errorobs, obsrec)) {
2354 0 : err += errorobs;
2355 : }
2356 21 : obslocation = mh.asMPosition();
2357 21 : }
2358 35 : err += verify(isSingleDish);
2359 35 : }
2360 0 : catch (AipsError &x) {
2361 0 : err = err + x.getMesg() + "\n";
2362 0 : }
2363 :
2364 35 : if (err.length() > 0) throw(AipsError("Invalid Image Parameter set : " + err));
2365 35 : }
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 35 : String SynthesisParamsImage::verify(const casacore::Bool isSingleDish) const
2401 : {
2402 35 : LogIO os;
2403 35 : String err;
2404 :
2405 35 : if(imageName=="") {err += "Please supply an image name\n";}
2406 :
2407 35 : if(imsize.nelements() != 2){err += "imsize must be a vector of 2 Ints\n";}
2408 35 : if(cellsize.nelements() != 2) {err += "cellsize must be a vector of 2 Quantities\n";}
2409 35 : 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 35 : if((mode=="mfs") && nchan > 1){
2418 0 : err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";
2419 : }
2420 :
2421 38 : if(! stokes.matches("I") && ! stokes.matches("Q") &&
2422 3 : ! stokes.matches("U") && ! stokes.matches("V") &&
2423 3 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2424 3 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2425 1 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2426 1 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2427 0 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2428 38 : ! 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 35 : 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 35 : Int imxnew = SynthesisUtilMethods::getOptimumSize(imsize[0]);
2468 35 : Int imynew = SynthesisUtilMethods::getOptimumSize(imsize[1]);
2469 35 : 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 70 : return err;
2483 35 : }// 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 145 : void SynthesisParamsImage::setDefaults()
2496 : {
2497 : // Image definition parameters
2498 145 : imageName = String("");
2499 145 : imsize.resize(2); imsize.set(100);
2500 145 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2501 145 : stokes="I";
2502 145 : phaseCenter=MDirection();
2503 145 : phaseCenterFieldId=-1;
2504 145 : projection=Projection::SIN;
2505 145 : useNCP=false;
2506 145 : startModel=Vector<String>(0);
2507 145 : freqFrameValid=True;
2508 145 : overwrite=false;
2509 : // PseudoI
2510 145 : pseudoi=false;
2511 :
2512 : // Spectral coordinates
2513 145 : nchan=1;
2514 145 : mode="mfs";
2515 145 : start="";
2516 145 : step="";
2517 145 : chanStart=0;
2518 145 : 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 145 : freqStart=Quantity(0,"");
2524 145 : freqStep=Quantity(0,"");
2525 145 : velStart=Quantity(0,"");
2526 145 : velStep=Quantity(0,"");
2527 145 : veltype=String("radio");
2528 145 : restFreq.resize(0);
2529 145 : refFreq = Quantity(0,"Hz");
2530 145 : frame = "";
2531 145 : freqFrame=MFrequency::LSRK;
2532 145 : sysvel="";
2533 145 : sysvelframe="";
2534 145 : sysvelvalue=Quantity(0.0,"m/s");
2535 145 : nTaylorTerms=1;
2536 145 : deconvolver="hogbom";
2537 : ///csysRecord=Record();
2538 145 : }
2539 :
2540 21 : Record SynthesisParamsImage::toRecord() const
2541 : {
2542 21 : Record impar;
2543 21 : impar.define("imagename", imageName);
2544 21 : impar.define("imsize", imsize);
2545 21 : Vector<String> cells(2);
2546 21 : cells[0] = QuantityToString( cellsize[0] );
2547 21 : cells[1] = QuantityToString( cellsize[1] );
2548 21 : impar.define("cell", cells );
2549 21 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2550 21 : else{impar.define("stokes", stokes);}
2551 21 : impar.define("nchan", nchan);
2552 21 : impar.define("nterms", nTaylorTerms);
2553 21 : impar.define("deconvolver",deconvolver);
2554 21 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2555 21 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2556 21 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2557 :
2558 21 : impar.define("specmode", mode );
2559 : // start and step can be one of these types
2560 21 : 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 21 : impar.define("start", start );
2578 : }
2579 21 : 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 21 : 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 21 : impar.define("veltype", veltype);
2606 21 : if (restfreqRecord.nfields() != 0 )
2607 : {
2608 0 : impar.defineRecord("restfreq", restfreqRecord);
2609 : }
2610 : else
2611 : {
2612 21 : Vector<String> rfs( restFreq.nelements() );
2613 21 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2614 21 : impar.define("restfreq", rfs);
2615 21 : }
2616 : //impar.define("reffreq", QuantityToString(refFreq));
2617 : //reffreq
2618 21 : if( reffreqRecord.nfields() != 0 )
2619 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2620 : else
2621 21 : { impar.define("reffreq",reffreq); }
2622 : //impar.define("reffreq", reffreq );
2623 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2624 21 : impar.define("outframe", frame );
2625 : //sysvel
2626 21 : if( sysvelRecord.nfields() != 0 )
2627 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2628 : else
2629 21 : { impar.define("sysvel", sysvel );}
2630 21 : impar.define("sysvelframe", sysvelframe );
2631 :
2632 21 : impar.define("restart",overwrite );
2633 21 : impar.define("freqframevalid", freqFrameValid);
2634 21 : impar.define("startmodel", startModel );
2635 :
2636 21 : if( csysRecord.isDefined("coordsys") )
2637 : {
2638 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2639 21 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2640 21 : impar.define("imshape", imshape);
2641 : }
2642 : // else cout << " NO CSYS INFO to write to output record " << endl;
2643 : ///Now save obslocation
2644 21 : Record tmprec;
2645 21 : String err;
2646 21 : MeasureHolder mh(obslocation);
2647 21 : if(mh.toRecord(err, tmprec)){
2648 21 : impar.defineRecord("obslocation_rec", tmprec);
2649 : }
2650 : else{
2651 0 : throw(AipsError("failed to save obslocation to record"));
2652 : }
2653 42 : return impar;
2654 21 : }
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 14 : 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 14 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2674 14 : vi2.originChunks();
2675 14 : vi2.origin();
2676 : /// This version uses the new vi2/vb2
2677 : // get the first ms for multiple MSes
2678 : //MeasurementSet msobj=vi2.ms();
2679 14 : Int fld=vb->fieldId()(0);
2680 :
2681 : //handling first ms only
2682 14 : Double gfreqmax=-1.0;
2683 14 : Double gdatafend=-1.0;
2684 14 : Double gdatafstart=1e14;
2685 14 : Double gfreqmin=1e14;
2686 14 : Vector<Int> spwids0;
2687 14 : Int j=0;
2688 14 : Int minfmsid=0;
2689 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2690 14 : Double imStartFreq=getCubeImageStartFreq();
2691 14 : std::vector<Int> sourceMsWithStartFreq;
2692 :
2693 :
2694 28 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2695 : //auto forMS0=chansel.find(0);
2696 14 : map<Int, Vector<Int> > spwsels=forMS0->second;
2697 14 : Int nspws=spwsels.size();
2698 14 : Vector<Int> spwids(nspws);
2699 14 : Vector<Int> nChannels(nspws);
2700 14 : Vector<Int> firstChannels(nspws);
2701 : //Vector<Int> channelIncrement(nspws);
2702 :
2703 14 : Int k=0;
2704 32 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2705 18 : spwids[k]=it->first;
2706 18 : nChannels[k]=(it->second)[0];
2707 18 : firstChannels[k]=(it->second)[1];
2708 : }
2709 14 : if(j==0) {
2710 14 : spwids0.resize();
2711 14 : 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 14 : Double freqmin=0, freqmax=0;
2723 14 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2724 :
2725 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2726 14 : 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 14 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2732 : //cerr << "after " << datafstart << " " << datafend << endl;
2733 14 : if((datafstart > datafend) || !status)
2734 0 : throw(AipsError("spw selection failed"));
2735 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2736 :
2737 14 : if (mode=="cubedata") {
2738 0 : freqmin = datafstart;
2739 0 : freqmax = datafend;
2740 : }
2741 14 : 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 28 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2766 28 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2767 : //cerr << "after " << freqmin << " " << freqmax << endl;
2768 : }
2769 :
2770 :
2771 :
2772 :
2773 14 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2774 14 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2775 14 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2776 14 : 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 14 : 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 14 : }
2802 14 : 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 14 : MeasurementSet msobj = *mss[minfmsid];
2809 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2810 28 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2811 14 : }
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 14 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2844 : MeasurementSet& msobj,
2845 : Vector<Int> spwids, Int fld,
2846 : Double freqmin, Double freqmax,
2847 : Double datafstart, Double datafend )
2848 : {
2849 28 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2850 :
2851 14 : CoordinateSystem csys;
2852 14 : 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 14 : MSColumns msc(msobj);
2865 14 : String telescop = msc.observation().telescopeName()(0);
2866 14 : MEpoch obsEpoch = msc.timeMeas()(0);
2867 14 : MPosition obsPosition;
2868 14 : 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 14 : MDirection phaseCenterToUse = phaseCenter;
2880 :
2881 14 : if( phaseCenterFieldId != -1 )
2882 : {
2883 3 : MSFieldColumns msfield(msobj.field());
2884 3 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2885 : {
2886 3 : if(trackSource){
2887 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2888 : }
2889 : else{
2890 3 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2891 : }
2892 : }
2893 : else
2894 : {
2895 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2896 : }
2897 3 : }
2898 : // Setup Phase center if it is specified only by field id.
2899 :
2900 : /////////////////// Direction Coordinates
2901 14 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2902 : // Normalize correctly
2903 14 : MVAngle ra=mvPhaseCenter.get()(0);
2904 14 : ra(0.0);
2905 :
2906 14 : MVAngle dec=mvPhaseCenter.get()(1);
2907 14 : Vector<Double> refCoord(2);
2908 14 : refCoord(0)=ra.get().getValue();
2909 14 : refCoord(1)=dec;
2910 14 : Vector<Double> refPixel(2);
2911 14 : refPixel(0) = Double(imsize[0]/2);
2912 14 : refPixel(1) = Double(imsize[1]/2);
2913 :
2914 14 : Vector<Double> deltas(2);
2915 14 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2916 14 : deltas(1)=cellsize[1].get("rad").getValue();
2917 14 : Matrix<Double> xform(2,2);
2918 14 : xform=0.0;xform.diagonal()=1.0;
2919 :
2920 14 : Vector<Double> projparams(2);
2921 14 : projparams = 0.0;
2922 14 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2923 14 : Projection projTo( projection.type(), projparams );
2924 :
2925 : DirectionCoordinate
2926 14 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2927 : // projection,
2928 : projTo,
2929 42 : refCoord(0), refCoord(1),
2930 42 : deltas(0), deltas(1),
2931 : xform,
2932 28 : 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 14 : obslocation=obsPosition;
2940 14 : ObsInfo myobsinfo;
2941 14 : myobsinfo.setTelescope(telescop);
2942 14 : myobsinfo.setPointingCenter(mvPhaseCenter);
2943 14 : myobsinfo.setObsDate(obsEpoch);
2944 14 : 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 14 : if(spwids.nelements()==0)
2960 : {
2961 0 : Int nspw=msc.spectralWindow().nrow();
2962 0 : spwids.resize(nspw);
2963 0 : indgen(spwids);
2964 : }
2965 14 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2966 14 : Vector<Double> dataChanFreq, dataChanWidth;
2967 14 : std::vector<std::vector<Int> > averageWhichChan;
2968 14 : std::vector<std::vector<Int> > averageWhichSPW;
2969 14 : std::vector<std::vector<Double> > averageChanFrac;
2970 :
2971 14 : if(spwids.nelements()==1)
2972 : {
2973 12 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2974 12 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
2975 : }
2976 : else
2977 : {
2978 2 : 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 14 : Double minDataFreq = min(dataChanFreq);
2985 14 : 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 14 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3037 14 : String cubemode;
3038 14 : if ( qrestfreq.getValue("Hz")==0 )
3039 : {
3040 14 : MSDopplerUtil msdoppler(msobj);
3041 14 : Vector<Double> restfreqvec;
3042 14 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3043 14 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3044 14 : 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 14 : }
3062 : Double refPix;
3063 14 : Vector<Double> chanFreq;
3064 14 : Vector<Double> chanFreqStep;
3065 14 : String specmode;
3066 :
3067 14 : 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 14 : 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 14 : Bool nonLinearFreq(false);
3083 14 : String veltype_p=veltype;
3084 14 : veltype_p.upcase();
3085 42 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3086 42 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3087 : {
3088 0 : nonLinearFreq= true;
3089 : }
3090 :
3091 14 : SpectralCoordinate mySpectral;
3092 : Double stepf;
3093 14 : if(!nonLinearFreq)
3094 : //TODO: velocity mode default start case (use last channels?)
3095 : {
3096 14 : Double startf=chanFreq[0];
3097 : //Double stepf=chanFreqStep[0];
3098 14 : if(chanFreq.nelements()==1)
3099 : {
3100 7 : stepf=chanFreqStep[0];
3101 : }
3102 : else
3103 : {
3104 7 : stepf=chanFreq[1]-chanFreq[0];
3105 : }
3106 14 : Double restf=qrestfreq.getValue("Hz");
3107 : //stepf=9e8;
3108 14 : 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 14 : 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 14 : 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 28 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3129 14 : 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 14 : uInt nrestfreq = restFreq.nelements();
3159 14 : 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 14 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3177 14 : if(whichStokes.nelements()==0)
3178 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3179 14 : StokesCoordinate myStokes(whichStokes);
3180 :
3181 : ////////////////// Build Full coordinate system.
3182 :
3183 : //CoordinateSystem csys;
3184 14 : csys.addCoordinate(myRaDec);
3185 14 : csys.addCoordinate(myStokes);
3186 14 : csys.addCoordinate(mySpectral);
3187 14 : csys.setObsInfo(myobsinfo);
3188 :
3189 : //store back csys to impars record
3190 : //cerr<<"save csys to csysRecord..."<<endl;
3191 14 : if(csysRecord.isDefined("coordsys"))
3192 0 : csysRecord.removeField("coordsys");
3193 14 : csys.save(csysRecord,"coordsys");
3194 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3195 : // imshape
3196 14 : imshape.resize(4);
3197 14 : imshape[0] = imsize[0];
3198 14 : imshape[1] = imsize[0];
3199 14 : imshape[2] = whichStokes.nelements();
3200 14 : 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 14 : } // end of else when coordsys record is not defined...
3234 :
3235 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3236 :
3237 28 : return csys;
3238 14 : }
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 14 : 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 14 : String inStart, inStep;
3531 14 : specmode = findSpecMode(mode);
3532 14 : String freqframe;
3533 14 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3534 28 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3535 :
3536 14 : refPix=0.0;
3537 14 : Bool descendingfreq(false);
3538 14 : Bool descendingoutfreq(false);
3539 :
3540 14 : if( mode.contains("cube") )
3541 : {
3542 13 : String restfreq=QuantityToString(qrestfreq);
3543 : // use frame from input start or width in MFreaquency or MRadialVelocity
3544 13 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3545 : // emit warning here if qmframe is used
3546 : //
3547 13 : inStart = start;
3548 13 : inStep = step;
3549 13 : if( specmode=="channel" )
3550 : {
3551 13 : inStart = String::toString(chanStart);
3552 13 : inStep = String::toString(chanStep);
3553 : // negative step -> descending channel indices
3554 13 : 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 13 : if (inStep=='0') inStep="";
3585 :
3586 13 : MRadialVelocity mSysVel;
3587 13 : Quantity qVel;
3588 : MRadialVelocity::Types mRef;
3589 13 : if(mode!="cubesource")
3590 : {
3591 :
3592 :
3593 13 : 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 13 : 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 13 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3626 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3627 26 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3628 13 : << LogIO::POST;
3629 13 : ostringstream ostr;
3630 13 : ostr << " phaseCenter='" << phaseCenter;
3631 13 : os << String(ostr)<<"' ";
3632 :
3633 : Double dummy; // dummy variable - weightScale is not used here
3634 26 : 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 13 : veltype,
3651 : verbose,
3652 : mSysVel
3653 : );
3654 :
3655 13 : if( nchan==-1 )
3656 : {
3657 7 : nchan = chanFreq.nelements();
3658 7 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3659 : }
3660 :
3661 13 : 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 13 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3668 26 : << LogIO::POST;
3669 :
3670 13 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3671 0 : descendingoutfreq = true;
3672 : }
3673 :
3674 : //if (descendingfreq && !descendingoutfreq) {
3675 26 : if ((specmode=="channel" && descendingfreq==1)
3676 26 : || (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 13 : }
3688 1 : else if ( mode=="mfs" ) {
3689 1 : chanFreq.resize(1);
3690 1 : chanFreqStep.resize(1);
3691 : //chanFreqStep[0] = freqmax - freqmin;
3692 1 : Double freqmean = (freqmin + freqmax)/2;
3693 1 : 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 1 : chanFreq[0] = refFreq.getValue("Hz");
3700 : // Set the new reffreq to be the refPix (CAS-9518)
3701 1 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3702 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3703 1 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3704 : }
3705 :
3706 1 : if( nchan==-1 ) nchan=1;
3707 1 : 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 14 : return true;
3719 :
3720 14 : }//getImFreq
3721 : /////////////////////////
3722 14 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3723 14 : Double inStartFreq=-1.0;
3724 14 : String checkspecmode("");
3725 14 : if(mode.contains("cube")) {
3726 13 : checkspecmode = findSpecMode(mode);
3727 : }
3728 14 : if(checkspecmode!="") {
3729 13 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3730 13 : if(checkspecmode=="channel") {
3731 13 : 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 14 : return inStartFreq;
3748 :
3749 14 : }
3750 :
3751 27 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3752 : {
3753 27 : String specmode;
3754 27 : specmode="channel";
3755 27 : 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 52 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3761 26 : !( velStep.getValue()==0.0)) {
3762 0 : specmode="velocity";
3763 : }
3764 52 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3765 26 : !(freqStep.getValue()==0.0)) {
3766 0 : specmode="frequency";
3767 : }
3768 : }
3769 27 : return specmode;
3770 0 : }
3771 :
3772 :
3773 42 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3774 : {
3775 42 : Vector<Int> whichStokes(0);
3776 60 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3777 27 : stokes=="RR" ||stokes=="LL" ||
3778 60 : stokes=="XX" || stokes=="YY" ) {
3779 39 : whichStokes.resize(1);
3780 39 : whichStokes(0)=Stokes::type(stokes);
3781 : }
3782 9 : else if(stokes=="IV" || stokes=="IQ" ||
3783 6 : stokes=="RRLL" || stokes=="XXYY" ||
3784 6 : stokes=="QU" || stokes=="UV"){
3785 3 : whichStokes.resize(2);
3786 :
3787 3 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3788 3 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3789 3 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3790 3 : 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 42 : return whichStokes;
3810 0 : }// decidenpolplanes
3811 :
3812 28 : IPosition SynthesisParamsImage::shp() const
3813 : {
3814 28 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3815 :
3816 28 : 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 56 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3822 : }
3823 :
3824 14 : Record SynthesisParamsImage::getcsys() const
3825 : {
3826 14 : 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 110 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3883 : {
3884 110 : setDefaults();
3885 110 : }
3886 :
3887 110 : SynthesisParamsGrid::~SynthesisParamsGrid()
3888 : {
3889 110 : }
3890 :
3891 :
3892 35 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3893 : {
3894 35 : setDefaults();
3895 :
3896 35 : String err("");
3897 :
3898 : try {
3899 35 : err += readVal( inrec, String("imagename"), imageName );
3900 :
3901 : // FTMachine parameters
3902 35 : err += readVal( inrec, String("gridder"), gridder );
3903 35 : err += readVal( inrec, String("padding"), padding );
3904 35 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3905 35 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3906 35 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3907 35 : err += readVal( inrec, String("convfunc"), convFunc );
3908 :
3909 35 : err += readVal( inrec, String("vptable"), vpTable );
3910 :
3911 :
3912 :
3913 : // convert 'gridder' to 'ftmachine' and 'mtype'
3914 35 : ftmachine = "gridft";
3915 35 : mType = "default";
3916 35 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
3917 13 : ftmachine = "gridft";
3918 : }
3919 22 : 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 22 : else if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes==1))
3925 0 : {ftmachine=="gridft";}
3926 :
3927 22 : else if (gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" ) {
3928 15 : ftmachine = "mosaicft";
3929 : }
3930 :
3931 7 : else if (gridder=="imagemosaic") {
3932 0 : mType = "imagemosaic";
3933 0 : if (wprojplanes>1 || wprojplanes==-1) {
3934 0 : ftmachine = "wprojectft";
3935 : }
3936 : }
3937 :
3938 7 : else if (gridder=="awproject" || gridder=="awprojectft" || gridder=="awp") {
3939 0 : ftmachine = "awprojectft";
3940 : }
3941 :
3942 7 : else if (gridder=="singledish") {
3943 :
3944 7 : ftmachine = "sd";
3945 : }
3946 : else{
3947 0 : ftmachine=gridder;
3948 0 : ftmachine.downcase();
3949 :
3950 : }
3951 :
3952 35 : String deconvolver;
3953 35 : err += readVal( inrec, String("deconvolver"), deconvolver );
3954 35 : if (deconvolver=="mtmfs") {
3955 1 : mType = "multiterm"; // Takes precedence over imagemosaic
3956 : }
3957 :
3958 : // facets
3959 35 : err += readVal( inrec, String("facets"), facets );
3960 : // chanchunks
3961 35 : err += readVal( inrec, String("chanchunks"), chanchunks );
3962 :
3963 : // Spectral interpolation
3964 35 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
3965 : // Track moving source ?
3966 35 : err += readVal( inrec, String("distance"), distance );
3967 35 : err += readVal( inrec, String("tracksource"), trackSource );
3968 35 : err += readVal( inrec, String("trackdir"), trackDir );
3969 :
3970 : // The extra params for WB-AWP
3971 35 : err += readVal( inrec, String("aterm"), aTermOn );
3972 35 : err += readVal( inrec, String("psterm"), psTermOn );
3973 35 : err += readVal( inrec, String("mterm"), mTermOn );
3974 35 : err += readVal( inrec, String("wbawp"), wbAWP );
3975 35 : err += readVal( inrec, String("cfcache"), cfCache );
3976 35 : err += readVal( inrec, String("usepointing"), usePointing );
3977 35 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3978 35 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3979 35 : err += readVal( inrec, String("conjbeams"), conjBeams );
3980 35 : err += readVal( inrec, String("computepastep"), computePAStep );
3981 35 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3982 :
3983 : // The extra params for single-dish
3984 35 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3985 35 : err += readVal( inrec, String("convertfirst"), convertFirst );
3986 35 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3987 35 : err += readVal( inrec, String("convsupport"), convSupport );
3988 35 : err += readVal( inrec, String("truncate"), truncateSize );
3989 35 : err += readVal( inrec, String("gwidth"), gwidth );
3990 35 : err += readVal( inrec, String("jwidth"), jwidth );
3991 35 : err += readVal( inrec, String("minweight"), minWeight );
3992 35 : 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 35 : if (ftmachine=="awprojectft" && cfCache=="") {
3998 0 : cfCache = imageName + ".cf";
3999 : }
4000 :
4001 35 : if ( ftmachine=="awprojectft" &&
4002 35 : 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 35 : err += verify();
4011 :
4012 35 : } catch(AipsError &x) {
4013 0 : err = err + x.getMesg() + "\n";
4014 0 : }
4015 :
4016 35 : if (err.length()>0) {
4017 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4018 : }
4019 :
4020 35 : }
4021 :
4022 35 : String SynthesisParamsGrid::verify() const
4023 : {
4024 35 : 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 35 : if ( imageName == "" ) {
4032 0 : err += "Please supply an image name\n";
4033 : }
4034 :
4035 44 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4036 71 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4037 64 : (ftmachine != "mawprojectft") &&
4038 7 : (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 70 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4049 35 : ( 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 35 : 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 35 : 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 35 : 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 35 : 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 35 : 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 35 : 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 35 : 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 35 : if ( ftmachine == "sd" ) {
4115 14 : if ( convertFirst != "always" and
4116 14 : 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 35 : return err;
4124 0 : }
4125 :
4126 145 : void SynthesisParamsGrid::setDefaults()
4127 : {
4128 145 : imageName="";
4129 : // FTMachine parameters
4130 : //ftmachine="GridFT";
4131 145 : ftmachine="gridft";
4132 145 : gridder=ftmachine;
4133 145 : padding=1.2;
4134 145 : useAutoCorr=false;
4135 145 : useDoublePrec=true;
4136 145 : wprojplanes=1;
4137 145 : convFunc="SF";
4138 145 : vpTable="";
4139 :
4140 : // facets
4141 145 : facets=1;
4142 :
4143 : // chanchunks
4144 145 : chanchunks=1;
4145 :
4146 : // Spectral Axis interpolation
4147 145 : interpolation=String("nearest");
4148 :
4149 : //mosaic use pointing
4150 145 : usePointing=false;
4151 : // Moving phase center ?
4152 145 : distance=Quantity(0,"m");
4153 145 : trackSource=false;
4154 145 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4155 :
4156 : // The extra params for WB-AWP
4157 145 : aTermOn = true;
4158 145 : psTermOn = true;
4159 145 : mTermOn = false;
4160 145 : wbAWP = true;
4161 145 : cfCache = "";
4162 145 : usePointing = false;
4163 145 : pointingOffsetSigDev.resize(0);
4164 : // pointingOffsetSigDev.set(30.0);
4165 145 : doPBCorr = true;
4166 145 : conjBeams = true;
4167 145 : computePAStep=360.0;
4168 145 : rotatePAStep=5.0;
4169 :
4170 : // extra params for single-dish
4171 145 : pointingDirCol = "";
4172 145 : convertFirst = "never";
4173 145 : skyPosThreshold = 0.0;
4174 145 : convSupport = -1;
4175 145 : truncateSize = Quantity(-1.0);
4176 145 : gwidth = Quantity(-1.0);
4177 145 : jwidth = Quantity(-1.0);
4178 145 : minWeight = 0.0;
4179 145 : clipMinMax = False;
4180 :
4181 : // Mapper type
4182 145 : mType = String("default");
4183 :
4184 145 : }
4185 :
4186 21 : Record SynthesisParamsGrid::toRecord() const
4187 : {
4188 21 : Record gridpar;
4189 :
4190 21 : gridpar.define("imagename", imageName);
4191 : // FTMachine params
4192 21 : gridpar.define("padding", padding);
4193 21 : gridpar.define("useautocorr",useAutoCorr );
4194 21 : gridpar.define("usedoubleprec", useDoublePrec);
4195 21 : gridpar.define("wprojplanes", wprojplanes);
4196 21 : gridpar.define("convfunc", convFunc);
4197 21 : gridpar.define("vptable", vpTable);
4198 :
4199 21 : gridpar.define("facets", facets);
4200 21 : gridpar.define("chanchunks", chanchunks);
4201 :
4202 21 : gridpar.define("interpolation",interpolation);
4203 :
4204 21 : gridpar.define("distance", QuantityToString(distance));
4205 21 : gridpar.define("tracksource", trackSource);
4206 21 : gridpar.define("trackdir", MDirectionToString( trackDir ));
4207 :
4208 21 : gridpar.define("aterm",aTermOn );
4209 21 : gridpar.define("psterm",psTermOn );
4210 21 : gridpar.define("mterm",mTermOn );
4211 21 : gridpar.define("wbawp", wbAWP);
4212 21 : gridpar.define("cfcache", cfCache);
4213 21 : gridpar.define("usepointing",usePointing );
4214 21 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
4215 21 : gridpar.define("dopbcorr", doPBCorr);
4216 21 : gridpar.define("conjbeams",conjBeams );
4217 21 : gridpar.define("computepastep", computePAStep);
4218 21 : gridpar.define("rotatepastep", rotatePAStep);
4219 :
4220 21 : gridpar.define("pointingcolumntouse", pointingDirCol );
4221 21 : gridpar.define("convertfirst", convertFirst );
4222 21 : gridpar.define("skyposthreshold", skyPosThreshold );
4223 21 : gridpar.define("convsupport", convSupport );
4224 21 : gridpar.define("truncate", QuantityToString(truncateSize) );
4225 21 : gridpar.define("gwidth", QuantityToString(gwidth) );
4226 21 : gridpar.define("jwidth", QuantityToString(jwidth) );
4227 21 : gridpar.define("minweight", minWeight );
4228 21 : gridpar.define("clipminmax", clipMinMax );
4229 :
4230 21 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
4231 : /// else gridpar.define("deconvolver","singleterm");
4232 :
4233 21 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
4234 21 : else gridpar.define("gridder", gridder);
4235 :
4236 : // gridpar.define("mtype", mType);
4237 :
4238 21 : return gridpar;
4239 0 : }
4240 :
4241 :
4242 :
4243 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
4244 :
4245 : /////////////////////// Deconvolver Parameters
4246 :
4247 32 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4248 : {
4249 32 : setDefaults();
4250 32 : }
4251 :
4252 32 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4253 : {
4254 32 : }
4255 :
4256 :
4257 31 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4258 : {
4259 31 : setDefaults();
4260 :
4261 31 : String err("");
4262 :
4263 : try
4264 : {
4265 :
4266 31 : err += readVal( inrec, String("imagename"), imageName );
4267 31 : 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 31 : if( inrec.isDefined("startmodel") )
4273 : {
4274 31 : if( inrec.dataType("startmodel")==TpString )
4275 : {
4276 5 : String onemodel;
4277 5 : err += readVal( inrec, String("startmodel"), onemodel );
4278 5 : if( onemodel.length()>0 )
4279 : {
4280 0 : startModel.resize(1);
4281 0 : startModel[0] = onemodel;
4282 : }
4283 5 : else {startModel.resize();}
4284 5 : }
4285 52 : else if( inrec.dataType("startmodel")==TpArrayString ||
4286 26 : inrec.dataType("startmodel")==TpArrayBool)
4287 : {
4288 26 : 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 31 : err += readVal( inrec, String("id"), deconvolverId );
4297 31 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4298 :
4299 31 : err += readVal( inrec, String("scales"), scales );
4300 31 : err += readVal( inrec, String("scalebias"), scalebias );
4301 :
4302 31 : err += readVal( inrec, String("usemask"), maskType );
4303 31 : if( maskType=="auto-thresh" )
4304 : {
4305 0 : autoMaskAlgorithm = "thresh";
4306 : }
4307 31 : else if( maskType=="auto-thesh2" )
4308 : {
4309 0 : autoMaskAlgorithm = "thresh2";
4310 : }
4311 31 : else if( maskType=="auto-onebox" )
4312 : {
4313 0 : autoMaskAlgorithm = "onebox";
4314 : }
4315 31 : else if( maskType=="user" || maskType=="pb" )
4316 : {
4317 31 : autoMaskAlgorithm = "";
4318 : }
4319 :
4320 :
4321 31 : if( inrec.isDefined("mask") )
4322 : {
4323 31 : if( inrec.dataType("mask")==TpString )
4324 : {
4325 17 : err+= readVal( inrec, String("mask"), maskString );
4326 : }
4327 14 : else if( inrec.dataType("mask")==TpArrayString )
4328 : {
4329 14 : err+= readVal( inrec, String("mask"), maskList );
4330 : }
4331 : }
4332 :
4333 31 : if( inrec.isDefined("pbmask") )
4334 : {
4335 31 : err += readVal( inrec, String("pbmask"), pbMask );
4336 : }
4337 31 : if( inrec.isDefined("maskthreshold") )
4338 : {
4339 31 : if( inrec.dataType("maskthreshold")==TpString )
4340 : {
4341 31 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4342 : //deal with the case a string is a float value without unit
4343 31 : Quantity testThresholdString;
4344 31 : Quantity::read(testThresholdString,maskThreshold);
4345 31 : if( testThresholdString.getUnit()=="" )
4346 : {
4347 31 : if(testThresholdString.getValue()<1.0)
4348 : {
4349 31 : fracOfPeak = testThresholdString.getValue();
4350 31 : maskThreshold=String("");
4351 : }
4352 : }
4353 31 : }
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 31 : if( inrec.isDefined("maskresolution") )
4372 : {
4373 31 : if( inrec.dataType("maskresolution")==TpString )
4374 : {
4375 31 : err += readVal(inrec, String("maskresolution"), maskResolution );
4376 : //deal with the case a string is a float value without unit
4377 31 : Quantity testResolutionString;
4378 31 : Quantity::read(testResolutionString,maskResolution);
4379 31 : if( testResolutionString.getUnit()=="" )
4380 : {
4381 31 : maskResByBeam = testResolutionString.getValue();
4382 31 : maskResolution=String("");
4383 : }
4384 31 : }
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 31 : if( inrec.isDefined("nmask") )
4398 : {
4399 31 : if( inrec.dataType("nmask")==TpInt )
4400 : {
4401 31 : err+= readVal(inrec, String("nmask"), nMask );
4402 : }
4403 : else
4404 : {
4405 0 : err+= "nmask must be an integer\n";
4406 : }
4407 : }
4408 31 : if( inrec.isDefined("autoadjust") )
4409 : {
4410 26 : if( inrec.dataType("autoadjust")==TpBool )
4411 : {
4412 26 : 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 31 : if (inrec.isDefined("fusedthreshold"))
4421 : {
4422 31 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4423 31 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4424 : else
4425 0 : err += "fusedthreshold must be a float or double";
4426 : }
4427 31 : if (inrec.isDefined("specmode"))
4428 : {
4429 31 : if(inrec.dataType("specmode") == TpString)
4430 31 : 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 31 : if (inrec.isDefined("largestscale"))
4436 : {
4437 31 : if (inrec.dataType("largestscale") == TpInt)
4438 31 : 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 31 : if( inrec.isDefined("sidelobethreshold"))
4444 : {
4445 31 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4446 : {
4447 31 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4448 : }
4449 : else
4450 : {
4451 0 : err+= "sidelobethreshold must be a float or double";
4452 : }
4453 : }
4454 :
4455 31 : if( inrec.isDefined("noisethreshold"))
4456 : {
4457 31 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4458 : {
4459 31 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4460 : }
4461 : else
4462 : {
4463 0 : err+= "noisethreshold must be a float or double";
4464 : }
4465 : }
4466 31 : if( inrec.isDefined("lownoisethreshold"))
4467 : {
4468 31 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4469 : {
4470 31 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4471 : }
4472 : else
4473 : {
4474 0 : err+= "lownoisethreshold must be a float or double";
4475 : }
4476 : }
4477 31 : if( inrec.isDefined("negativethreshold"))
4478 : {
4479 31 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4480 : {
4481 31 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4482 : }
4483 : else
4484 : {
4485 0 : err+= "negativethreshold must be a float or double";
4486 : }
4487 : }
4488 31 : if( inrec.isDefined("smoothfactor"))
4489 : {
4490 31 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4491 : {
4492 31 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4493 : }
4494 : else
4495 : {
4496 0 : err+= "smoothfactor must be a float or double";
4497 : }
4498 : }
4499 31 : if( inrec.isDefined("minbeamfrac"))
4500 : {
4501 31 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4502 : {
4503 31 : 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 31 : if( inrec.isDefined("cutthreshold"))
4517 : {
4518 31 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4519 : {
4520 31 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4521 : }
4522 : else {
4523 0 : err+= "cutthreshold must be a float or double";
4524 : }
4525 : }
4526 31 : if( inrec.isDefined("growiterations"))
4527 : {
4528 31 : if (inrec.dataType("growiterations")==TpInt) {
4529 31 : err+= readVal(inrec, String("growiterations"), growIterations );
4530 : }
4531 : else {
4532 0 : err+= "growiterations must be an integer\n";
4533 : }
4534 : }
4535 31 : if( inrec.isDefined("dogrowprune"))
4536 : {
4537 31 : if (inrec.dataType("dogrowprune")==TpBool) {
4538 31 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4539 : }
4540 : else {
4541 0 : err+= "dogrowprune must be a bool\n";
4542 : }
4543 : }
4544 31 : if( inrec.isDefined("minpercentchange"))
4545 : {
4546 31 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4547 31 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4548 : }
4549 : else {
4550 0 : err+= "minpercentchange must be a float or double";
4551 : }
4552 : }
4553 31 : if( inrec.isDefined("verbose"))
4554 : {
4555 31 : if (inrec.dataType("verbose")==TpBool ) {
4556 31 : err+= readVal(inrec, String("verbose"), verbose);
4557 : }
4558 : else {
4559 0 : err+= "verbose must be a bool";
4560 : }
4561 : }
4562 31 : if( inrec.isDefined("fastnoise"))
4563 : {
4564 31 : if (inrec.dataType("fastnoise")==TpBool ) {
4565 31 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4566 : }
4567 : else {
4568 0 : err+= "fastnoise must be a bool";
4569 : }
4570 : }
4571 31 : if( inrec.isDefined("nsigma") )
4572 : {
4573 31 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4574 31 : 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 31 : if( inrec.isDefined("noRequireSumwt") )
4587 : {
4588 26 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4589 26 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4590 : }
4591 : else {
4592 0 : err+= "noRequireSumwt must be a bool";
4593 : }
4594 : }
4595 31 : if( inrec.isDefined("fullsummary") )
4596 : {
4597 31 : if (inrec.dataType("fullsummary")==TpBool) {
4598 31 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4599 : }
4600 : else {
4601 0 : err+= "fullsummary must be a bool";
4602 : }
4603 : }
4604 31 : if( inrec.isDefined("restoringbeam") )
4605 : {
4606 5 : String errinfo("");
4607 : try {
4608 :
4609 5 : 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 5 : 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 5 : }// if isdefined(restoringbeam)
4658 :
4659 31 : if( inrec.isDefined("interactive") )
4660 : {
4661 31 : if( inrec.dataType("interactive")==TpBool )
4662 31 : {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 31 : err += verify();
4670 :
4671 : }
4672 0 : catch(AipsError &x)
4673 : {
4674 0 : err = err + x.getMesg() + "\n";
4675 0 : }
4676 :
4677 31 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4678 :
4679 31 : }
4680 :
4681 31 : String SynthesisParamsDeconv::verify() const
4682 : {
4683 31 : String err;
4684 :
4685 31 : 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 31 : if( maskString.length()>0 )
4689 : {
4690 4 : File fp( imageName+".mask" );
4691 4 : 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 4 : }
4693 :
4694 31 : if( pbMask >= 1.0)
4695 0 : {err += "pbmask must be < 1.0 \n"; }
4696 31 : else if( pbMask < 0.0)
4697 0 : {err += "pbmask must be a positive value \n"; }
4698 :
4699 31 : 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 31 : if ( fracOfPeak >= 1.0)
4708 0 : {err += "fracofpeak must be < 1.0 \n"; }
4709 31 : else if ( fracOfPeak < 0.0)
4710 0 : {err += "fracofpeak must be a positive value \n"; }
4711 :
4712 31 : return err;
4713 0 : }
4714 :
4715 63 : void SynthesisParamsDeconv::setDefaults()
4716 : {
4717 63 : imageName="";
4718 63 : algorithm="hogbom";
4719 63 : startModel=Vector<String>(0);
4720 63 : deconvolverId=0;
4721 63 : nTaylorTerms=1;
4722 63 : scales.resize(1); scales[0]=0.0;
4723 63 : scalebias=0.6;
4724 63 : maskType="none";
4725 63 : maskString="";
4726 63 : maskList.resize(1); maskList[0]="";
4727 63 : pbMask=0.0;
4728 63 : autoMaskAlgorithm="thresh";
4729 63 : maskThreshold="";
4730 63 : maskResolution="";
4731 63 : fracOfPeak=0.0;
4732 63 : nMask=0;
4733 63 : interactive=false;
4734 63 : autoAdjust=False;
4735 63 : fusedThreshold = 0.0;
4736 63 : specmode="mfs";
4737 63 : largestscale = -1;
4738 63 : }
4739 :
4740 26 : Record SynthesisParamsDeconv::toRecord() const
4741 : {
4742 26 : Record decpar;
4743 :
4744 26 : decpar.define("imagename", imageName);
4745 26 : decpar.define("deconvolver", algorithm);
4746 26 : decpar.define("startmodel",startModel);
4747 26 : decpar.define("id",deconvolverId);
4748 26 : decpar.define("nterms",nTaylorTerms);
4749 26 : decpar.define("scales",scales);
4750 26 : decpar.define("scalebias",scalebias);
4751 26 : decpar.define("usemask",maskType);
4752 26 : decpar.define("fusedthreshold", fusedThreshold);
4753 26 : decpar.define("specmode", specmode);
4754 26 : decpar.define("largestscale", largestscale);
4755 26 : if( maskList.nelements()==1 && maskList[0]=="")
4756 : {
4757 12 : decpar.define("mask",maskString);
4758 : }
4759 : else {
4760 14 : decpar.define("mask",maskList);
4761 : }
4762 26 : decpar.define("pbmask",pbMask);
4763 26 : if (fracOfPeak > 0.0)
4764 : {
4765 0 : decpar.define("maskthreshold",fracOfPeak);
4766 : }
4767 : else
4768 : {
4769 26 : decpar.define("maskthreshold",maskThreshold);
4770 : }
4771 26 : decpar.define("maskresolution",maskResolution);
4772 26 : decpar.define("nmask",nMask);
4773 26 : decpar.define("autoadjust",autoAdjust);
4774 26 : decpar.define("sidelobethreshold",sidelobeThreshold);
4775 26 : decpar.define("noisethreshold",noiseThreshold);
4776 26 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4777 26 : decpar.define("negativethreshold",negativeThreshold);
4778 26 : decpar.define("smoothfactor",smoothFactor);
4779 26 : decpar.define("minbeamfrac",minBeamFrac);
4780 26 : decpar.define("cutthreshold",cutThreshold);
4781 26 : decpar.define("growiterations",growIterations);
4782 26 : decpar.define("dogrowprune",doGrowPrune);
4783 26 : decpar.define("minpercentchange",minPercentChange);
4784 26 : decpar.define("verbose", verbose);
4785 26 : decpar.define("fastnoise", fastnoise);
4786 26 : decpar.define("interactive",interactive);
4787 26 : decpar.define("nsigma",nsigma);
4788 26 : decpar.define("noRequireSumwt",noRequireSumwt);
4789 26 : decpar.define("fullsummary",fullsummary);
4790 :
4791 26 : return decpar;
4792 0 : }
4793 :
4794 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4795 :
4796 :
4797 : } //# NAMESPACE CASA - END
4798 :
|
