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 37 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 37 : }
89 :
90 37 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 37 : }
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 6120 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 6120 : Int N=vb.nRows(),M=-1;
108 6120 : for(Int i=0;i<N;i++)
109 : {
110 6120 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 6120 : {M++;break;}
112 : }
113 6120 : 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 34 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 34 : Int n=npix;
121 :
122 34 : if( n%2 !=0 ){ n+= 1; }
123 :
124 34 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 272 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 238 : if( fac[k]>5 )
129 : {
130 0 : val = fac[k];
131 0 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 0 : fac[k] = val;
133 : }
134 : }
135 34 : newlarge=product(fac);
136 34 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 0 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 34 : return newlarge;
141 34 : }
142 :
143 : // Return the prime factors of the given number
144 34 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 34 : Vector<uInt> factors;
147 :
148 34 : Int lastresult = n;
149 34 : Int sqlast = int(sqrt(n))+1;
150 :
151 34 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 34 : Int c=2;
153 : while(1)
154 : {
155 272 : if( lastresult == 1 || c > sqlast ) { break; }
156 238 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 238 : if( c>sqlast){ c=lastresult; break; }
160 238 : if( lastresult % c == 0 ) { break; }
161 0 : c += 1;
162 : }
163 238 : factors.resize( factors.nelements()+1, true );
164 238 : factors[ factors.nelements()-1 ] = c;
165 238 : lastresult /= c;
166 : }
167 34 : 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 34 : return factors;
185 0 : }
186 :
187 :
188 0 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 0 : if (psfcutoff >=1.0 || psfcutoff<=0.0)
193 : {
194 0 : os << "psfcutoff must be >0 and <1" << LogIO::WARN;
195 0 : return false;
196 : }
197 :
198 0 : std::shared_ptr<SIImageStore> imstore;
199 0 : if( nterms>1 )
200 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 0 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 0 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 0 : imstore->printBeamSet();
210 :
211 0 : imstore->releaseLocks();
212 :
213 0 : return true;
214 0 : }
215 :
216 :
217 :
218 : // Evaluate a polynomial from Taylor coefficients and expand to a Cube.
219 : // Since this funcion is primarily for use with mtmfs_via_cube, this step applies only to
220 : // the multiterm.model.ttx images and the cube.model image cube.
221 0 : Bool SynthesisUtilMethods::taylorCoeffsToCube(const String& cubename,const String& mtname, const Int nterms, const String& reffreq)
222 : {
223 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "taylorCoeffsToCube"));
224 :
225 : // Set up imstores
226 0 : CountedPtr<SIImageStore> cube_imstore;
227 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
228 :
229 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
230 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
231 :
232 : // Check that .model exists.
233 : try{
234 0 : cube_imstore->model();
235 0 : for(Int i=0;i<nterms;i++)
236 0 : mt_imstore->model(i);
237 : }
238 0 : catch(AipsError &x)
239 : {
240 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\nModel images must exist on disk." ));
241 0 : }
242 :
243 : // Get/check shapes
244 0 : IPosition cube_shp( cube_imstore->model()->shape() );
245 0 : IPosition mt_shp( mt_imstore->model(0)->shape() );
246 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
247 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
248 : }
249 :
250 : // Read reference frequency
251 0 : Quantity reffreq_qa;
252 0 : Quantity::read( reffreq_qa, reffreq );
253 0 : Double refval = reffreq_qa.getValue("Hz");
254 :
255 : //cout << "ref freq : " << refval << endl;
256 :
257 : //Get the frequency list for the cube
258 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
259 0 : Vector<Double> freqlist( cube_shp[3] );
260 :
261 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
262 : {
263 0 : if( csys.type(i) == Coordinate::SPECTRAL )
264 : {
265 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
266 :
267 0 : for(Int ch=0;ch<cube_shp[3];ch++)
268 : {
269 : Double freq;
270 0 : Bool ret = speccoord.toWorld( freq, ch );
271 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
272 0 : freqlist[ch] = freq;
273 :
274 : //cout << "freq " << ch << " is " << freq << endl;
275 : }
276 :
277 0 : }
278 : }
279 :
280 :
281 : // Reset the Cube values to zero.
282 : //cube_imstore->model()->set(0.0);
283 :
284 : //For each pol, do the Taylor-to-Cube calculation.
285 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
286 : {
287 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(nterms);
288 0 : for(Int i=0;i<nterms;i++)
289 : {
290 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
291 0 : 0, cube_shp[3],
292 0 : pol, cube_shp[2],
293 0 : *mt_imstore->model(i) );
294 : }
295 :
296 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
297 : {
298 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
299 0 : chan, cube_shp[3],
300 0 : pol, cube_shp[2],
301 0 : *cube_imstore->model() );
302 :
303 0 : Double wt = (freqlist[chan] - refval) / refval;
304 :
305 0 : cube_subim->set(0.0);
306 0 : for(Int tt=0;tt<nterms;tt++)
307 : {
308 0 : Double fac = pow(wt,tt);
309 0 : LatticeExpr<Float> oneterm = LatticeExpr<Float>( *cube_subim + (fac) * (*mt_subims[tt])) ;
310 0 : cube_subim->copyData(oneterm);
311 0 : }
312 :
313 :
314 0 : }//for chan
315 0 : }// for pol
316 :
317 0 : return True;
318 :
319 0 : }//end of func
320 :
321 :
322 : // Calculate the RHS of the Normal equations for a linear least squares fit of a Taylor polynomial (per pixel).
323 : // This function is primarily for use with mtmfs_via_cube, and may be used for the PSF (2nterms-1), the residual (nterms)
324 : // and the primary beam (nterms=1).
325 : // imtype=0 : PSF with 2nterms-1 terms
326 : // imtype=1 : residual with nterms terms
327 : // imtype=2 : pb with 1 term
328 : // imtype=3 : sumwt with 2nterms-1 terms
329 0 : Bool SynthesisUtilMethods::cubeToTaylorSum(const String& cubename,const String& mtname, const Int nterms, const String& reffreq, const Int imtype, const Float pblimit)
330 : {
331 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "cubeToTaylorSum"));
332 :
333 : //cout << "imtype : " << imtype << endl;
334 0 : if(imtype <0 || imtype >3)
335 : {
336 0 : throw( AipsError("cubeToTaylorSum currently only supports 'psf','residual','pb', 'sumwt' options"));
337 : }
338 :
339 : // Set up imstores
340 0 : CountedPtr<SIImageStore> cube_imstore;
341 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
342 :
343 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
344 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
345 : //For psf avg pb has to be done already
346 : // doing residual too for sensitivity this is independent of beam spectral index removal
347 0 : Float maxPB=1.0;
348 0 : if(imtype < 2){
349 0 : LatticeExprNode elnod( max( *(mt_imstore->pb(0)) ) );
350 0 : maxPB=elnod.getFloat();
351 0 : if(maxPB == 0.0){
352 0 : throw(AipsError("Programmers error: should do tt psf images after making average PB"));
353 :
354 : }
355 :
356 0 : }
357 : // cerr << "imtype " << imtype << " MAX PB " << maxPB << endl;
358 : // If dopsf=True, calculate 2n-1 terms.
359 0 : Int out_nterms=nterms; // for residual
360 0 : if(imtype==0 || imtype==3){out_nterms=2 * nterms - 1;} // the psfs fill the upper triangle of the Hessian with 2 nterms-1 elements. Also sumwt.
361 0 : if(imtype==2){out_nterms=1;} // For the PB, for mtmfs_via_cube, we need only tt0. Later, if we need all terms to calculate PB alpha, then change this to nterms, and add the invHesian math (elsewhere) to later convert the RHS vector into the coefficients.
362 :
363 0 : CountedPtr <ImageInterface<Float> > use_cube, use_mt;
364 : // If dopsf=True, check that .psf cube and mt's exist. If dopsf=False, check residual images.
365 : try{
366 0 : switch(imtype)
367 : {
368 0 : case 0: use_cube=cube_imstore->psf();break;
369 0 : case 1: use_cube=cube_imstore->residual();break;
370 0 : case 2: use_cube=cube_imstore->pb();break;
371 0 : case 3: use_cube=cube_imstore->sumwt();break;
372 : }
373 0 : cube_imstore->sumwt();
374 0 : for(Int i=0;i<out_nterms;i++)
375 : {
376 0 : switch(imtype)
377 : {
378 0 : case 0:mt_imstore->psf(i);break;
379 0 : case 1:mt_imstore->residual(i);break;
380 0 : case 2:mt_imstore->pb(i);break;
381 0 : case 3:mt_imstore->sumwt(i);break;
382 : }
383 : }
384 : }
385 0 : catch(AipsError &x)
386 : {
387 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n " + imtype + " images must exist on disk." ));
388 0 : }
389 :
390 : // Get/check shapes ( Assume that the PSF always exists in the imstore... A valid assumption in the context of mtmfs_via_cube )
391 0 : IPosition cube_shp( cube_imstore->psf()->shape() );
392 0 : IPosition mt_shp( mt_imstore->psf(0)->shape() );
393 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
394 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
395 : }
396 :
397 : // Read reference frequency
398 0 : Quantity reffreq_qa;
399 0 : Quantity::read( reffreq_qa, reffreq );
400 0 : Double refval = reffreq_qa.getValue("Hz");
401 :
402 : //cout << "ref freq : " << refval << endl;
403 :
404 : //Get the frequency list for the cube
405 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
406 0 : Vector<Double> freqlist( cube_shp[3] );
407 :
408 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
409 : {
410 0 : if( csys.type(i) == Coordinate::SPECTRAL )
411 : {
412 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
413 :
414 0 : for(Int ch=0;ch<cube_shp[3];ch++)
415 : {
416 : Double freq;
417 0 : Bool ret = speccoord.toWorld( freq, ch );
418 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
419 0 : freqlist[ch] = freq;
420 : // cout << "freq " << ch << " is " << freq << endl;
421 : }
422 :
423 0 : }
424 : }
425 :
426 :
427 : // Reset the Taylor Sum values to zero.
428 0 : for(Int i=0;i<out_nterms;i++)
429 : {
430 0 : switch(imtype)
431 : {
432 0 : case 0:mt_imstore->psf(i)->set(0.0);break;
433 0 : case 1:mt_imstore->residual(i)->set(0.0);break;
434 0 : case 2:mt_imstore->pb(i)->set(0.0);break;
435 0 : case 3:mt_imstore->sumwt(i)->set(0.0);break;
436 : }
437 : }
438 :
439 : // Get the sumwt spectrum.
440 0 : Array<Float> lsumwt;
441 0 : cube_imstore->sumwt()->get(lsumwt, False);
442 :
443 : // Sum the weights ( or just use accumulate...)
444 0 : LatticeExprNode msum( sum( *cube_imstore->sumwt() ) );
445 0 : Float wtsum = msum.getFloat();
446 :
447 : //cerr << "perchansumwt : shape "<< lsumwt.shape() << " " << lsumwt << " sumwt "<< wtsum << endl;
448 :
449 : //Float wtsum = cube_shp[3]; // This is sum of weights, if all weights are 1.0
450 :
451 : //For each pol, do the Cube-To-Taylor calculation.
452 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
453 : {
454 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(out_nterms);
455 0 : for(Int i=0;i<out_nterms;i++)
456 : {
457 0 : switch(imtype)
458 : {
459 0 : case 0:use_mt=mt_imstore->psf(i);break;
460 0 : case 1:use_mt=mt_imstore->residual(i);break;
461 0 : case 2:use_mt=mt_imstore->pb(i);break;
462 0 : case 3:use_mt=mt_imstore->sumwt(i);break;
463 : }
464 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
465 0 : 0, cube_shp[3],
466 0 : pol, cube_shp[2],
467 0 : *use_mt );
468 : }
469 :
470 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
471 : {
472 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
473 0 : chan, cube_shp[3],
474 0 : pol, cube_shp[2],
475 0 : *use_cube);
476 0 : if(imtype < 2){
477 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
478 0 : chan, cube_shp[3],
479 0 : pol, cube_shp[2],
480 0 : *(cube_imstore->pb()));
481 0 : CountedPtr<ImageInterface<Float> > tmplat = new TempImage<Float>(cube_subim->shape(), cube_subim->coordinates());
482 0 : tmplat->copyData(LatticeExpr<Float>((*pb_subim) *(*cube_subim)));
483 0 : cube_subim = tmplat;
484 :
485 0 : }
486 :
487 0 : IPosition pos(4,0,0,pol,chan);
488 :
489 0 : Double wt = (freqlist[chan] - refval) / refval;
490 :
491 0 : for(Int tt=0;tt<out_nterms;tt++)
492 : {
493 0 : Double fac = pow(wt,tt);
494 : //cerr << "BEF accum " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
495 0 : LatticeExpr<Float> eachterm = LatticeExpr<Float>( (*mt_subims[tt]) + ((fac) * (*cube_subim) * lsumwt(pos))) ;
496 0 : mt_subims[tt]->copyData(eachterm);
497 : //cerr <<" AFT accum : chan " << chan << " tt " << tt << " fac " << fac << " lsumwt " << lsumwt(pos) << " pos " << pos << " max " << max(mt_subims[tt]->get()) << endl;
498 0 : }
499 :
500 :
501 0 : }//for chan
502 :
503 :
504 : // Divide by sum of weights.
505 0 : for(Int tt=0;tt<out_nterms;tt++)
506 : {
507 : //cerr << "bef div : tt " << tt << " : " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
508 :
509 0 : LatticeExpr<Float> eachterm;
510 0 : if (imtype < 2) {
511 0 : eachterm = LatticeExpr<Float>( iif( (*(mt_imstore->pb(0))) > pblimit , (*mt_subims[tt]) / wtsum/(*(mt_imstore->pb(0))), 0.0));
512 : }
513 : else{
514 0 : eachterm = LatticeExpr<Float>( (*mt_subims[tt]) / wtsum ) ;
515 : }
516 0 : mt_subims[tt]->copyData(eachterm);
517 :
518 0 : mt_subims[tt]->flush();
519 : // cerr << "aft div : " << max(mt_subims[tt]->get()) << endl;
520 0 : }
521 :
522 0 : }// for pol
523 :
524 :
525 : // Set the T/F mask, for PB images. Without this, the PB is fully masked, for aproj /mosaic gridders.
526 0 : if( imtype==2 )
527 : {
528 0 : mt_imstore->removeMask( mt_imstore->pb(0) );
529 : {
530 : //MSK//
531 0 : LatticeExpr<Bool> pbmask( iif( *mt_imstore->pb(0) > fabs(pblimit) , True , False ) );
532 : //MSK//
533 0 : mt_imstore->createMask( pbmask, mt_imstore->pb(0) );
534 0 : mt_imstore->pb(0)->pixelMask().unlock();
535 0 : }
536 :
537 : }
538 :
539 0 : return True;
540 :
541 0 : }//end of func
542 :
543 :
544 0 : Bool SynthesisUtilMethods::removeFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
545 : {
546 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "removeFreqDepPB"));
547 :
548 : // Set up imstores
549 0 : CountedPtr<SIImageStore> cube_imstore;
550 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
551 :
552 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
553 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
554 :
555 : try{
556 0 : cube_imstore->residual(); // Residual Cube
557 0 : cube_imstore->pb(); // PB cube
558 0 : mt_imstore->pb(0); // avgPB in the tt0 pb.
559 : }
560 0 : catch(AipsError &x)
561 : {
562 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n Residual cube, PB cube, and multiterm PB.tt0 must exist on disk." ));
563 0 : }
564 :
565 : // Get/check shapes
566 0 : IPosition cube_shp( cube_imstore->residual()->shape() );
567 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
568 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
569 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
570 : }
571 :
572 : //For each pol, do the freq-dep PB math.//////////////////////////
573 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
574 : {
575 :
576 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
577 0 : 0, cube_shp[3],
578 0 : pol, cube_shp[2],
579 0 : (*mt_imstore->pb(0)) );
580 :
581 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
582 0 : Float mtpbmaxval = mtpbmax.getFloat();
583 0 : if(mtpbmaxval <=0.0){os << LogIO::WARN << "pb.tt0 max is < or = zero" << LogIO::POST;}
584 :
585 :
586 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
587 : {
588 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
589 0 : chan, cube_shp[3],
590 0 : pol, cube_shp[2],
591 0 : *cube_imstore->residual() );
592 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
593 0 : chan, cube_shp[3],
594 0 : pol, cube_shp[2],
595 0 : *cube_imstore->pb() );
596 :
597 0 : LatticeExprNode pbmax( max( *pb_subim ) );
598 0 : Float pbmaxval = pbmax.getFloat();
599 0 : if( pbmaxval<=0.0 )
600 : {
601 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
602 : }
603 : else
604 : {
605 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*mt_subim)*(*(cube_subim))/(*(pb_subim)) , 0.0 ) );
606 0 : cube_subim->copyData( thepbcor );
607 0 : }// if not zero
608 :
609 0 : }//for chan
610 0 : }// for pol
611 :
612 0 : return True;
613 :
614 0 : }//end of func
615 :
616 :
617 :
618 0 : Bool SynthesisUtilMethods::applyFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
619 : {
620 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "applyFreqDepPB"));
621 :
622 : // Set up imstores
623 0 : CountedPtr<SIImageStore> cube_imstore;
624 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
625 :
626 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
627 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
628 :
629 : try{
630 0 : cube_imstore->model(); // Model Cube
631 0 : cube_imstore->pb(); // PB cube
632 0 : mt_imstore->pb(0); // avgPB in the tt0 pb.
633 : }
634 0 : catch(AipsError &x)
635 : {
636 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n Model cube, PB cube, and multiterm PB.tt0 must exist on disk." ));
637 0 : }
638 :
639 : // Get/check shapes
640 0 : IPosition cube_shp( cube_imstore->model()->shape() );
641 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
642 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
643 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
644 : }
645 :
646 : //For each pol, do the freq-dep PB math.//////////////////////////
647 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
648 : {
649 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
650 0 : 0, cube_shp[3],
651 0 : pol, cube_shp[2],
652 0 : (*mt_imstore->pb(0)) );
653 :
654 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
655 0 : Float mtpbmaxval = mtpbmax.getFloat();
656 0 : if(mtpbmaxval <=0.0)
657 0 : {os << LogIO::SEVERE << "pb.tt0 max is < or = zero. Cannot divide model image ! ERROR" << LogIO::POST;}
658 : else
659 : {
660 :
661 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
662 : {
663 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
664 0 : chan, cube_shp[3],
665 0 : pol, cube_shp[2],
666 0 : *cube_imstore->model() );
667 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
668 0 : chan, cube_shp[3],
669 0 : pol, cube_shp[2],
670 0 : *cube_imstore->pb() );
671 :
672 :
673 0 : LatticeExprNode pbmax( max( *pb_subim ) );
674 0 : Float pbmaxval = pbmax.getFloat();
675 0 : if( pbmaxval<=0.0 )
676 : {
677 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
678 : }
679 : else
680 : {
681 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*(cube_subim)) *(*(pb_subim)) / (*mt_subim) , 0.0 ) );
682 0 : cube_subim->copyData( thepbcor );
683 0 : }// if not zero
684 :
685 0 : }//for chan
686 : }// if mtpb >0
687 0 : }// for pol
688 :
689 0 : return True;
690 :
691 0 : }//end of func
692 :
693 :
694 :
695 :
696 : /***make a record of synthesisimager::weight parameters***/
697 17 : 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 17 : Record outRec;
705 17 : outRec.define("type", type);
706 17 : outRec.define("rmode", rmode);
707 17 : Record quantRec;
708 17 : QuantumHolder(noise).toRecord(quantRec);
709 17 : outRec.defineRecord("noise", quantRec);
710 17 : outRec.define("robust", robust);
711 17 : QuantumHolder(fieldofview).toRecord(quantRec);
712 17 : outRec.defineRecord("fieldofview", quantRec);
713 17 : outRec.define("npixels", npixels);
714 17 : outRec.define("multifield", multiField);
715 17 : outRec.define("usecubebriggs", useCubeBriggs);
716 17 : outRec.define("filtertype", filtertype);
717 17 : QuantumHolder(filterbmaj).toRecord(quantRec);
718 17 : outRec.defineRecord("filterbmaj", quantRec);
719 17 : QuantumHolder(filterbmin).toRecord(quantRec);
720 17 : outRec.defineRecord("filterbmin", quantRec);
721 17 : QuantumHolder(filterbpa).toRecord(quantRec);
722 17 : outRec.defineRecord("filterbpa", quantRec);
723 17 : outRec.define("fracBW", fracBW);
724 :
725 34 : return outRec;
726 17 : }
727 0 : void SynthesisUtilMethods::getFromWeightRecord(String& type, String& rmode,
728 : Quantity& noise, Double& robust,
729 : Quantity& fieldofview,
730 : Int& npixels, Bool& multiField, Bool& useCubeBriggs,
731 : String& filtertype, Quantity& filterbmaj,
732 : Quantity& filterbmin, Quantity& filterbpa, Double& fracBW, const Record& inRec){
733 0 : QuantumHolder qh;
734 0 : String err;
735 0 : if(!inRec.isDefined("type"))
736 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
737 0 : inRec.get("type", type);
738 0 : inRec.get("rmode", rmode);
739 0 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
740 0 : throw(AipsError("Error in reading noise param"));
741 0 : noise=qh.asQuantity();
742 0 : inRec.get("robust", robust);
743 0 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
744 0 : throw(AipsError("Error in reading fieldofview param"));
745 0 : fieldofview=qh.asQuantity();
746 0 : inRec.get("npixels", npixels);
747 0 : inRec.get("multifield", multiField);
748 0 : inRec.get("usecubebriggs", useCubeBriggs);
749 0 : inRec.get("filtertype", filtertype);
750 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
751 0 : throw(AipsError("Error in reading filterbmaj param"));
752 0 : filterbmaj=qh.asQuantity();
753 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
754 0 : throw(AipsError("Error in reading filterbmin param"));
755 0 : filterbmin=qh.asQuantity();
756 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
757 0 : throw(AipsError("Error in reading filterbpa param"));
758 0 : filterbpa=qh.asQuantity();
759 0 : inRec.get("fracBW", fracBW);
760 :
761 :
762 :
763 0 : }
764 :
765 :
766 : /**
767 : * Get values from lines of a /proc/self/status file. For example:
768 : * 'VmRSS: 600 kB'
769 : * @param str line from status file
770 : * @return integer value (memory amount, etc.)
771 : */
772 0 : Int SynthesisUtilMethods::parseProcStatusLine(const std::string &str) {
773 0 : istringstream is(str);
774 0 : std::string token;
775 0 : is >> token;
776 0 : is >> token;
777 0 : Int value = stoi(token);
778 0 : return value;
779 0 : }
780 :
781 : /**
782 : * Produces a name for a 'memprofile' output file. For example:
783 : * casa.synthesis.imager.memprofile.23514.pc22555.hq.eso.org.20171209_120446.txt
784 : * (where 23514 is the PID passed as input parameter).
785 : *
786 : * @param pid PID of the process running the imager
787 : *
788 : * @return a longish 'memprofile' filename including PID, machine, timestamp, etc.
789 : **/
790 0 : String SynthesisUtilMethods::makeResourceFilename(int pid)
791 : {
792 0 : if (g_hostname.empty() or g_startTimestamp.empty()) {
793 : // TODO: not using HOST_NAME_MAX because of issues with __APPLE__
794 : // somehow tests tAWPFTM, tSynthesisImager, and tSynthesisImagerVi2 fail.
795 0 : const int strMax = 255;
796 : char hostname[strMax];
797 0 : gethostname(hostname, strMax);
798 0 : g_hostname = hostname;
799 :
800 0 : auto time = std::time(nullptr);
801 0 : auto gmt = std::gmtime(&time);
802 0 : const char* format = "%Y%m%d_%H%M%S";
803 : char timestr[strMax];
804 0 : std::strftime(timestr, strMax, format, gmt);
805 0 : g_startTimestamp = timestr;
806 : }
807 :
808 0 : return String("casa.synthesis.imager.memprofile." + String::toString(pid) +
809 0 : "." + g_hostname + "." + g_startTimestamp + ".txt");
810 : }
811 :
812 0 : Bool SynthesisUtilMethods::adviseChanSel(Double& freqStart, Double& freqEnd,
813 : const Double& freqStep, const MFrequency::Types& freqframe,
814 : Vector<Int>& spw, Vector<Int>& start,
815 : Vector<Int>& nchan, const String& ms, const String& ephemtab, const Int field_id, const Bool getFreqRange, const String spwselection){
816 :
817 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "adviseChanSel"));
818 0 : if(ms==String("")){
819 0 : throw(AipsError("Need a valid MS"));
820 : }
821 0 : spw.resize();
822 0 : start.resize();
823 0 : nchan.resize();
824 : try {
825 0 : if(!getFreqRange){
826 0 : Vector<Int> bnchan;
827 0 : Vector<Int> bstart;
828 0 : Vector<Int> bspw;
829 : Double fS, fE;
830 0 : fS=freqStart;
831 0 : fE=freqEnd;
832 0 : if(freqEnd < freqStart){
833 0 : fS=freqEnd;
834 0 : fE=freqStart;
835 : }
836 :
837 :
838 : {
839 :
840 0 : MeasurementSet elms(String(ms), TableLock(TableLock::AutoNoReadLocking), Table::Old);
841 0 : if(ephemtab != "" && freqframe == MFrequency::REST ){
842 0 : MSUtil::getSpwInSourceFreqRange(bspw, bstart, bnchan, elms, fS, fE, fabs(freqStep), ephemtab, field_id);
843 : }
844 : else
845 0 : MSUtil::getSpwInFreqRange(bspw, bstart, bnchan, elms, fS, fE, fabs(freqStep), freqframe, field_id);
846 0 : elms.relinquishAutoLocks(true);
847 :
848 0 : }
849 0 : spw=Vector<Int> (bspw);
850 0 : start=Vector<Int> (bstart);
851 0 : nchan=Vector<Int> (bnchan);
852 0 : }
853 : else{
854 :
855 : {
856 0 : MeasurementSet elms(ms, TableLock(TableLock::AutoNoReadLocking), Table::Old);
857 0 : MSSelection thisSelection;
858 0 : String spsel=spwselection;
859 0 : if(spsel=="")spsel="*";
860 0 : thisSelection.setSpwExpr(spsel);
861 0 : TableExprNode exprNode=thisSelection.toTableExprNode(&elms);
862 0 : Matrix<Int> chanlist=thisSelection.getChanList();
863 0 : if(chanlist.ncolumn() <3){
864 0 : freqStart=-1.0;
865 0 : freqEnd=-1.0;
866 0 : return false;
867 : }
868 0 : Vector<Int> elspw=chanlist.column(0);
869 0 : Vector<Int> elstart=chanlist.column(1);
870 0 : Vector<Int> elnchan=Vector<Int> (chanlist.column(2)-elstart)+1;
871 0 : if(ephemtab != "" ){
872 0 : const MSColumns mscol(ms);
873 0 : MEpoch ep=mscol.timeMeas()(0);
874 0 : Quantity sysvel;
875 0 : String ephemTable("");
876 0 : MDirection::Types mtype=MDirection::APP;
877 0 : MDirection mdir(mtype);
878 0 : if(Table::isReadable(ephemtab)){
879 0 : ephemTable=ephemtab;
880 : }
881 0 : else if(ephemtab=="TRACKFIELD"){
882 0 : ephemTable=(mscol.field()).ephemPath(field_id);
883 : }
884 0 : else if(MDirection::getType(mtype, ephemtab)){
885 0 : mdir=MDirection(mtype);
886 : }
887 :
888 0 : MSUtil::getFreqRangeAndRefFreqShift(freqStart, freqEnd, sysvel, ep, elspw, elstart, elnchan, elms, ephemTable , mdir, True);
889 :
890 0 : }
891 : else
892 0 : MSUtil::getFreqRangeInSpw(freqStart, freqEnd, elspw, elstart, elnchan, elms, freqframe, field_id);
893 0 : }
894 :
895 : }
896 :
897 :
898 :
899 :
900 0 : } catch (AipsError x) {
901 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
902 0 : << LogIO::POST;
903 0 : return false;
904 0 : }
905 0 : catch (...){
906 : os << LogIO::SEVERE << "Unknown exception handled"
907 0 : << LogIO::POST;
908 0 : return false;
909 :
910 0 : }
911 :
912 0 : return true;
913 :
914 0 : }
915 :
916 289 : 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 289 : Bool isOn = false;
924 289 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
925 289 : if (!isOn)
926 289 : return;
927 :
928 : // TODO: reorganize, use struct or something to hold and pass info over. ifdef lnx
929 0 : LogIO casalog( LogOrigin("SynthesisUtilMethods", "getResource", WHERE) );
930 :
931 : // To hold memory stats, in MB
932 0 : int vmRSS = -1, vmHWM = -1, vmSize = -1, vmPeak = -1, vmSwap = -1;
933 0 : pid_t pid = -1;
934 0 : int fdSize = -1;
935 :
936 : // TODO: this won't probably work on anything but linux
937 0 : ifstream procFile("/proc/self/status");
938 0 : if (procFile.is_open()) {
939 0 : std::string line;
940 0 : while (not procFile.eof()) {
941 0 : getline(procFile, line);
942 0 : const std::string startVmRSS = "VmRSS:";
943 0 : const std::string startVmWHM = "VmHWM:";
944 0 : const std::string startVmSize = "VmSize:";
945 0 : const std::string startVmPeak = "VmPeak:";
946 0 : const std::string startVmSwap = "VmSwap:";
947 0 : const std::string startFDSize = "FDSize:";
948 0 : const double KB_TO_MB = 1024.0;
949 0 : if (startVmRSS == line.substr(0, startVmRSS.size())) {
950 0 : vmRSS = parseProcStatusLine(line.c_str()) / KB_TO_MB;
951 0 : } else if (startVmWHM == line.substr(0, startVmWHM.size())) {
952 0 : vmHWM = parseProcStatusLine(line.c_str()) / KB_TO_MB;
953 0 : } else if (startVmSize == line.substr(0, startVmSize.size())) {
954 0 : vmSize = parseProcStatusLine(line.c_str()) / KB_TO_MB;
955 0 : } else if (startVmPeak == line.substr(0, startVmPeak.size())) {
956 0 : vmPeak = parseProcStatusLine(line.c_str()) / KB_TO_MB;
957 0 : } else if (startVmSwap == line.substr(0, startVmSwap.size())) {
958 0 : vmSwap = parseProcStatusLine(line.c_str()) / KB_TO_MB;
959 0 : } else if (startFDSize == line.substr(0, startFDSize.size())) {
960 0 : fdSize = parseProcStatusLine(line.c_str());
961 : }
962 0 : }
963 0 : procFile.close();
964 0 : }
965 :
966 0 : pid = getpid();
967 :
968 : struct rusage usage;
969 : struct timeval now;
970 0 : getrusage(RUSAGE_SELF, &usage);
971 0 : now = usage.ru_utime;
972 :
973 : // TODO: check if this works as expected when /proc/self/status is not there
974 : // Not clear at all if VmHWM and .ru_maxrss measure the same thing
975 : // Some alternative is needed for the other fields as well: VmSize, VMHWM, FDSize.
976 0 : if (vmHWM < 0) {
977 0 : vmHWM = usage.ru_maxrss;
978 : }
979 :
980 0 : ostringstream oss;
981 0 : oss << "PID: " << pid ;
982 0 : oss << " MemRSS (VmRSS): " << vmRSS << " MB.";
983 0 : oss << " VmWHM: " << vmHWM << " MB.";
984 0 : oss << " VirtMem (VmSize): " << vmSize << " MB.";
985 0 : oss << " VmPeak: " << vmPeak << " MB.";
986 0 : oss << " VmSwap: " << vmSwap << " MB.";
987 0 : oss << " ProcTime: " << now.tv_sec << '.' << now.tv_usec;
988 0 : oss << " FDSize: " << fdSize;
989 0 : oss << " [" << label << "] ";
990 0 : casalog << oss.str() << LogIO::NORMAL3 << LogIO::POST;
991 :
992 : // Write this to a file too...
993 : try {
994 0 : if (fname.empty()) {
995 0 : fname = makeResourceFilename(pid);
996 : }
997 0 : ofstream ofile(fname, ios::app);
998 0 : if (ofile.is_open()) {
999 0 : if (0 == ofile.tellp()) {
1000 : casalog << g_enableOptMemProfile << " is enabled, initializing output file for "
1001 : "imager profiling information (memory and run time): " << fname <<
1002 0 : LogIO::NORMAL << LogIO::POST;
1003 0 : ostringstream header;
1004 : header << "# PID, MemRSS_(VmRSS)_MB, VmWHM_MB, VirtMem_(VmSize)_MB, VmPeak_MB, "
1005 0 : "VmSwap_MB, ProcTime_sec, FDSize, label_checkpoint";
1006 0 : ofile << header.str() << '\n';
1007 0 : }
1008 0 : ostringstream line;
1009 0 : line << pid << ',' << vmRSS << ',' << vmHWM << ',' << vmSize << ','
1010 0 : << vmPeak << ','<< vmSwap << ',' << now.tv_sec << '.' << now.tv_usec << ','
1011 0 : << fdSize << ',' << '[' << label << ']';
1012 0 : ofile << line.str() << '\n';
1013 0 : ofile.close();
1014 0 : }
1015 0 : } catch(std::runtime_error &exc) {
1016 : casalog << "Could not write imager memory+runtime information into output file: "
1017 0 : << fname << LogIO::WARN << LogIO::POST;
1018 0 : }
1019 0 : }
1020 :
1021 :
1022 : // Data partitioning rules for CONTINUUM imaging
1023 : //
1024 : // ALL members of the selection parameters in selpars are strings
1025 : // ONLY. This methods reads the selection parameters from selpars
1026 : // and returns a partitioned Record with npart data selection
1027 : // entries.
1028 : //
1029 : // The algorithm used to do the partitioning along the TIME axis is
1030 : // as follows:
1031 : //
1032 : // for each MS in selpars
1033 : // - get the selection parameters
1034 : // - generate a selected MS
1035 : // - get number of rows in the selected MS
1036 : // - divide the rows in nparts
1037 : // - for each part
1038 : // - get compute rowBeg and rowEnd
1039 : // - modify rowEnd such that rowEnd points to the end of
1040 : // full integration data. This is done as follows:
1041 : // tRef = TIME(rowEnd);
1042 : // reduce rowEnd till TIME(rowEnd) != tRef
1043 : // - Construct a T0~T1 string
1044 : // - Fill it in the timeSelPerPart[msID][PartNo] array
1045 : //
1046 0 : Record SynthesisUtilMethods::continuumDataPartition(Record &selpars, const Int npart)
1047 : {
1048 0 : LogIO os( LogOrigin("SynthesisUtilMethods","continuumDataPartition",WHERE) );
1049 :
1050 0 : Record onepart, allparts;
1051 0 : Vector<Vector<String> > timeSelPerPart;
1052 0 : timeSelPerPart.resize(selpars.nfields());
1053 :
1054 : // Duplicate the entire input record npart times, with a specific partition id.
1055 : // Modify each sub-record according to partition id.
1056 0 : for (uInt msID=0;msID<selpars.nfields();msID++)
1057 : {
1058 0 : Record thisMS= selpars.subRecord(RecordFieldId("ms"+String::toString(msID)));
1059 0 : String msName = thisMS.asString("msname");
1060 0 : timeSelPerPart[msID].resize(npart,true);
1061 : //
1062 : // Make a selected MS and extract the time-column information
1063 : //
1064 0 : MeasurementSet ms(msName,TableLock(TableLock::AutoNoReadLocking), Table::Old),
1065 0 : selectedMS(ms);
1066 0 : MSInterface msI(ms); MSSelection msSelObj;
1067 0 : msSelObj.reset(msI,MSSelection::PARSE_NOW,
1068 : thisMS.asString("timestr"),
1069 : thisMS.asString("antenna"),
1070 : thisMS.asString("field"),
1071 : thisMS.asString("spw"),
1072 : thisMS.asString("uvdist"),
1073 : thisMS.asString("taql"),
1074 : "",//thisMS.asString("poln"),
1075 : thisMS.asString("scan"),
1076 : "",//thisMS.asString("array"),
1077 : thisMS.asString("state"),
1078 : thisMS.asString("obs")//observation
1079 : );
1080 0 : msSelObj.getSelectedMS(selectedMS);
1081 :
1082 : //--------------------------------------------------------------------
1083 : // Use the selectedMS to generate time selection strings per part
1084 : //
1085 : // Double Tint;
1086 0 : MSMainColumns mainCols(selectedMS);
1087 0 : Vector<rownr_t> rowNumbers = selectedMS.rowNumbers();
1088 0 : Int nRows=selectedMS.nrow(),
1089 0 : dRows=nRows/npart;
1090 0 : Int rowBegID=0, rowEndID=nRows-1;
1091 0 : Int rowBeg=rowNumbers[rowBegID], rowEnd=rowNumbers[rowEndID];
1092 : //cerr << "NRows, dRows, npart = " << nRows << " " << dRows << " " << npart << " " << rowBeg << " " << rowEnd << endl;
1093 :
1094 0 : rowEndID = rowBegID + dRows;
1095 :
1096 :
1097 0 : MVTime mvInt=mainCols.intervalQuant()(0);
1098 : //Time intT(mvInt.getTime());
1099 : // Tint = intT.modifiedJulianDay();
1100 :
1101 0 : Int partNo=0;
1102 : // The +1 in rowBeg and rowEnd below is because it appears
1103 : // that TaQL by default counts from 1, not 0.
1104 0 : while(rowEndID < nRows)
1105 : {
1106 : // rowBeg=rowNumbers[rowBegID]; rowEnd = rowNumbers[rowEndID];
1107 0 : rowBeg=rowBegID+1; rowEnd = rowEndID+1;
1108 0 : stringstream taql;
1109 0 : taql << "ROWNUMBER() >= " << rowBeg << " && ROWNUMBER() <= " << rowEnd;
1110 0 : timeSelPerPart[msID][partNo] = taql.str();
1111 :
1112 0 : if (partNo == npart - 1) break;
1113 0 : partNo++;
1114 0 : rowBegID = rowEndID+1;
1115 0 : rowEndID = min(rowBegID + dRows, nRows-1);
1116 0 : if (rowEndID == nRows-1) break;
1117 0 : }
1118 :
1119 : //rowBeg=rowNumbers[rowBegID]; rowEnd = rowNumbers[nRows-1];
1120 0 : stringstream taql;
1121 0 : rowBeg=rowBegID+1; rowEnd = nRows-1+1;
1122 0 : taql << "ROWNUMBER() >= " << rowBeg << " && ROWNUMBER() <= " << rowEnd;
1123 0 : timeSelPerPart[msID][partNo] = taql.str();
1124 : os << endl << "Rows = " << rowBeg << " " << rowEnd << " "
1125 0 : << "[P][M]: " << msID << ":" << partNo << " " << timeSelPerPart[msID][partNo]
1126 0 : << LogIO::POST;
1127 0 : }
1128 : //
1129 : // The time selection strings for all parts of the current
1130 : // msID are in timeSelPerPart.
1131 : //--------------------------------------------------------------------
1132 : //
1133 : // Now reverse the order of parts and ME loops. Create a Record
1134 : // per part, get the MS for thisPart. Put the associated time
1135 : // selection string in it. Add the thisMS to thisPart Record.
1136 : // Finally add thisPart Record to the allparts Record.
1137 : //
1138 0 : for(Int iPart=0; iPart<npart; iPart++)
1139 : {
1140 0 : Record thisPart;
1141 0 : thisPart.assign(selpars);
1142 0 : for (uInt msID=0; msID<selpars.nfields(); msID++)
1143 : {
1144 0 : Record thisMS = thisPart.subRecord(RecordFieldId("ms"+String::toString(msID)));
1145 :
1146 0 : thisMS.define("taql",timeSelPerPart[msID][iPart]);
1147 0 : thisPart.defineRecord(RecordFieldId("ms"+String::toString(msID)) , thisMS);
1148 0 : }
1149 0 : allparts.defineRecord(RecordFieldId(String::toString(iPart)), thisPart);
1150 0 : }
1151 : // cerr << allparts << endl;
1152 0 : return allparts;
1153 :
1154 : // for( Int part=0; part < npart; part++)
1155 : // {
1156 :
1157 : // onepart.assign(selpars);
1158 :
1159 :
1160 : // //-------------------------------------------------
1161 : // // WARNING : This is special-case code for two specific datasets
1162 : // for ( uInt msid=0; msid<selpars.nfields(); msid++)
1163 : // {
1164 : // Record onems = onepart.subRecord( RecordFieldId("ms"+String::toString(msid)) );
1165 : // String msname = onems.asString("msname");
1166 : // String spw = onems.asString("spw");
1167 : // if(msname.matches("DataTest/twopoints_twochan.ms"))
1168 : // {
1169 : // onems.define("spw", spw+":"+String::toString(part));
1170 : // }
1171 : // if(msname.matches("DataTest/point_twospws.ms"))
1172 : // {
1173 : // onems.define("spw", spw+":"+ (((Bool)part)?("10~19"):("0~9")) );
1174 : // }
1175 : // if(msname.matches("DataTest/reg_mawproject.ms"))
1176 : // {
1177 : // onems.define("scan", (((Bool)part)?("1~17"):("18~35")) );
1178 : // }
1179 : // onepart.defineRecord( RecordFieldId("ms"+String::toString(msid)) , onems );
1180 : // }// end ms loop
1181 : // //-------------------------------------------------
1182 :
1183 : // allparts.defineRecord( RecordFieldId(String::toString(part)), onepart );
1184 :
1185 : // }// end partition loop
1186 :
1187 : // return allparts;
1188 0 : }
1189 :
1190 :
1191 : // Data partitioning rules for CUBE imaging
1192 0 : Record SynthesisUtilMethods::cubeDataPartition(const Record &selpars, const Int npart,
1193 : const Double freqBeg, const Double freqEnd, const MFrequency::Types eltype)
1194 : {
1195 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataPartition",WHERE) );
1196 : // Temporary special-case code. Please replace with actual rules.
1197 0 : Vector<Double> fstart(npart);
1198 0 : Vector<Double> fend(npart);
1199 0 : Double step=(freqEnd-freqBeg)/Double(npart);
1200 0 : fstart(0)=freqBeg;
1201 0 : fend(0)=freqBeg+step;
1202 0 : for (Int k=1; k < npart; ++k){
1203 0 : fstart(k)=fstart(k-1)+step;
1204 0 : fend(k)=fend(k-1)+step;
1205 : }
1206 0 : return cubeDataPartition( selpars, fstart, fend, eltype );
1207 :
1208 0 : }
1209 :
1210 :
1211 0 : Record SynthesisUtilMethods::cubeDataImagePartition(const Record & selpars, const CoordinateSystem&
1212 : incsys, const Int npart, const Int nchannel,
1213 : Vector<CoordinateSystem>& outCsys,
1214 : Vector<Int>& outnChan){
1215 :
1216 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataImagePartition",WHERE) );
1217 0 : outnChan.resize(npart);
1218 0 : outCsys.resize(npart);
1219 0 : Int nomnchan=nchannel/npart;
1220 0 : outnChan.set(nomnchan);
1221 0 : nomnchan=nchannel%npart;
1222 0 : for (Int k=0; k < nomnchan; ++k)
1223 0 : outnChan[k]+=1;
1224 0 : Vector<Int> shp(0);
1225 : //shp(0)=20; shp(1)=20; shp(2)=1; shp(3)=outnChan[0];
1226 0 : Vector<Float> shift(4, 0.0);
1227 0 : Vector<Float> fac(4, 1.0);
1228 0 : Vector<Double> freqEnd(npart);
1229 0 : Vector<Double> freqStart(npart);
1230 0 : Float chanshift=0.0;
1231 0 : for (Int k =0; k <npart; ++k){
1232 0 : shift(3)=chanshift;
1233 : //cerr << k << " shift " << shift << endl;
1234 0 : outCsys[k]=incsys.subImage(shift, fac, shp);
1235 0 : freqStart[k]=SpectralImageUtil::worldFreq(outCsys[k], 0.0);
1236 0 : freqEnd[k]=SpectralImageUtil::worldFreq(outCsys[k], Double(outnChan[k]-1));
1237 0 : if(freqStart[k] > freqEnd[k]){
1238 0 : Double tmp=freqEnd[k];
1239 0 : freqEnd[k]=freqStart[k];
1240 0 : freqStart[k]=tmp;
1241 : }
1242 0 : chanshift+=Float(outnChan[k]);
1243 : }
1244 0 : MFrequency::Types eltype=incsys.spectralCoordinate(incsys.findCoordinate(Coordinate::SPECTRAL)).frequencySystem(true);
1245 :
1246 : //os << "freqStart="<<freqStart<<" freqend="<<freqEnd<< "eltype="<<eltype<<LogIO::POST;
1247 0 : Record rec=cubeDataPartition(selpars, freqStart, freqEnd, eltype);
1248 0 : for (Int k=0; k < npart ; ++k){
1249 0 : outCsys[k].save(rec.asrwRecord(String::toString(k)), "coordsys");
1250 0 : rec.asrwRecord(String::toString(k)).define("nchan", outnChan[k]);
1251 : }
1252 0 : return rec;
1253 0 : }
1254 :
1255 0 : Record SynthesisUtilMethods::cubeDataPartition(const Record &selpars, const Vector<Double>& freqBeg, const Vector<Double>&freqEnd, const MFrequency::Types eltype){
1256 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeDataPartition",WHERE) );
1257 0 : Record retRec;
1258 0 : Int npart=freqBeg.shape()(0);
1259 0 : for (Int k=0; k < npart; ++k){
1260 0 : Int nms=selpars.nfields();
1261 0 : Record partRec;
1262 0 : for(Int j=0; j < nms; ++j){
1263 0 : if(selpars.isDefined(String("ms"+String::toString(j)))){
1264 0 : Record msRec=selpars.asRecord(String("ms"+String::toString(j)));
1265 0 : if(!msRec.isDefined("msname"))
1266 0 : throw(AipsError("No msname key in ms record"));
1267 0 : String msname=msRec.asString("msname");
1268 0 : String userspw=msRec.isDefined("spw")? msRec.asString("spw") : "*";
1269 0 : String userfield=msRec.isDefined("field") ? msRec.asString("field") : "*";
1270 0 : String userstate=msRec.isDefined("state") ? msRec.asString("state") : "*";
1271 :
1272 0 : MeasurementSet elms(msname);
1273 0 : Record laSelection=elms.msseltoindex(userspw, userfield);
1274 0 : if (userfield=="") userfield="*";
1275 0 : MSSelection mssel;
1276 0 : mssel.setSpwExpr(userspw);
1277 0 : mssel.setFieldExpr(userfield);
1278 0 : mssel.setStateExpr(userstate);
1279 0 : TableExprNode exprNode = mssel.toTableExprNode(&elms);
1280 0 : Matrix<Int> spwsel=mssel.getChanList();
1281 0 : Vector<Int> fieldsel=mssel.getFieldList();
1282 : // case for scan intent specified
1283 0 : if (userstate!="*") {
1284 0 : MeasurementSet elselms((elms)(exprNode), &elms);
1285 0 : MSColumns tmpmsc(elselms);
1286 0 : Vector<Int> fldidv=tmpmsc.fieldId().getColumn();
1287 0 : if (fldidv.nelements()==0)
1288 0 : throw(AipsError("No field ids were selected, please check input parameters"));
1289 0 : std::set<Int> ufldids(fldidv.begin(),fldidv.end());
1290 0 : std::vector<Int> tmpv(ufldids.begin(), ufldids.end());
1291 0 : fieldsel.resize(tmpv.size());
1292 0 : uInt count=0;
1293 0 : for (std::vector<int>::const_iterator it=tmpv.begin();it != tmpv.end(); it++)
1294 : {
1295 0 : fieldsel(count) = *it;
1296 0 : count++;
1297 : }
1298 0 : }
1299 : //Matrix<Int> spwsel=laSelection.asArrayInt("channel");
1300 : //Vector<Int> fieldsel=laSelection.asArrayInt("field");
1301 0 : Vector<Int> freqSpw;
1302 0 : Vector<Int> freqStart;
1303 0 : Vector<Int> freqNchan;
1304 0 : String newspw;
1305 : try{
1306 0 : MSUtil::getSpwInFreqRange(freqSpw, freqStart, freqNchan, elms, freqBeg(k), freqEnd(k),0.0, eltype, fieldsel[0]);
1307 0 : newspw=mergeSpwSel(freqSpw, freqStart, freqNchan, spwsel);
1308 : //cerr << "try " << freqSpw << " " << freqStart << " " << freqNchan << endl;
1309 : }
1310 0 : catch(...){
1311 : //cerr << "In catch " << endl;
1312 0 : newspw="";
1313 0 : }
1314 : //String newspw=mergeSpwSel(freqSpw, freqStart, freqNchan, spwsel);
1315 0 : if(newspw=="") newspw="-1";
1316 0 : msRec.define("spw", newspw);
1317 0 : partRec.defineRecord(String("ms"+String::toString(j)),msRec);
1318 0 : }
1319 :
1320 : }
1321 0 : retRec.defineRecord(String::toString(k), partRec);
1322 0 : }
1323 :
1324 :
1325 :
1326 :
1327 0 : return retRec;
1328 0 : }
1329 :
1330 :
1331 0 : String SynthesisUtilMethods::mergeSpwSel(const Vector<Int>& fspw, const Vector<Int>& fstart, const Vector<Int>& fnchan, const Matrix<Int>& spwsel)
1332 : {
1333 0 : String retval="";
1334 : Int cstart, cend;
1335 0 : for(Int k=0; k < fspw.shape()(0); ++k){
1336 0 : cstart=fstart(k);
1337 0 : cend=fstart(k)+fnchan(k)-1;
1338 0 : for (Int j=0; j < spwsel.shape()(0); ++j){
1339 : //need to put this here as multiple rows can have the same spw
1340 0 : cstart=fstart(k);
1341 0 : cend=fstart(k)+fnchan(k)-1;
1342 0 : if(spwsel(j,0)==fspw[k]){
1343 0 : if(cstart < spwsel(j,1)) cstart=spwsel(j,1);
1344 0 : if(cend > spwsel(j,2)) cend= spwsel(j,2);
1345 0 : if(!retval.empty()) retval=retval+(",");
1346 0 : retval=retval+String::toString(fspw[k])+":"+String::toString(cstart)+"~"+String::toString(cend);
1347 : }
1348 : }
1349 : }
1350 :
1351 :
1352 :
1353 0 : return retval;
1354 0 : }
1355 :
1356 : // Image cube partitioning rules for CUBE imaging
1357 0 : Record SynthesisUtilMethods::cubeImagePartition(Record &impars, Int npart)
1358 : {
1359 0 : LogIO os( LogOrigin("SynthesisUtilMethods","cubeImagePartition",WHERE) );
1360 :
1361 0 : Record onepart, allparts;
1362 :
1363 : // Duplicate the entire input record npart times, with a specific partition id.
1364 : // Modify each sub-record according to partition id.
1365 0 : for( Int part=0; part < npart; part++)
1366 : {
1367 :
1368 0 : onepart.assign(impars);
1369 :
1370 : //-------------------------------------------------
1371 : // WARNING : This is special-case code for two specific datasets
1372 0 : for ( uInt imfld=0; imfld<impars.nfields(); imfld++)
1373 : {
1374 0 : Record onefld = onepart.subRecord( RecordFieldId(String::toString(imfld)) );
1375 0 : Int nchan = onefld.asInt("nchan");
1376 : //String freqstart = onems.asString("freqstart");
1377 :
1378 0 : onefld.define("nchan", nchan/npart);
1379 0 : onefld.define("freqstart", (((Bool)part)?("1.2GHz"):("1.0GHz")) );
1380 :
1381 0 : String imname = onefld.asString("imagename");
1382 0 : onefld.define("imagename", imname+".n"+String::toString(part));
1383 :
1384 0 : onepart.defineRecord( RecordFieldId( String::toString(imfld) ), onefld );
1385 0 : }// end ms loop
1386 : //-------------------------------------------------
1387 :
1388 0 : allparts.defineRecord( RecordFieldId(String::toString(part)), onepart );
1389 :
1390 : }// end partition loop
1391 :
1392 0 : return allparts;
1393 :
1394 :
1395 0 : }
1396 :
1397 0 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
1398 : {
1399 0 : MVAngle mvRa=direction.getAngle().getValue()(0);
1400 0 : MVAngle mvDec=direction.getAngle().getValue()(1);
1401 0 : ostringstream oos;
1402 0 : oos << " ";
1403 0 : Int widthRA=20;
1404 0 : Int widthDec=20;
1405 0 : oos.setf(ios::left, ios::adjustfield);
1406 0 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
1407 0 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
1408 0 : oos << " "
1409 0 : << MDirection::showType(direction.getRefPtr()->getType());
1410 0 : return String(oos);
1411 0 : }
1412 :
1413 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1414 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1415 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
1416 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1417 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1418 :
1419 : // Read String from Record
1420 1105 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1421 : {
1422 1105 : if( rec.isDefined( id ) )
1423 : {
1424 1020 : String inval("");
1425 1020 : if( rec.dataType( id )==TpString )
1426 1020 : { 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 1020 : if(inval.length()>0){val=inval;}
1432 1020 : return String("");
1433 : }
1434 0 : else { return String(id + " must be a string\n"); }
1435 1020 : }
1436 85 : else { return String("");}
1437 : }
1438 :
1439 : // Read Integer from Record
1440 272 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1441 : {
1442 272 : if( rec.isDefined( id ) )
1443 : {
1444 272 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1445 0 : else { return String(id + " must be an integer\n"); }
1446 : }
1447 0 : else { return String(""); }
1448 : }
1449 :
1450 : // Read Float from Record
1451 493 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1452 : {
1453 493 : if( rec.isDefined( id ) )
1454 : {
1455 459 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1456 459 : { rec.get( id , val ); return String(""); }
1457 0 : else { return String(id + " must be a float\n"); }
1458 : }
1459 34 : else { return String(""); }
1460 : }
1461 :
1462 : // Read Bool from Record
1463 527 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1464 : {
1465 527 : if( rec.isDefined( id ) )
1466 : {
1467 408 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1468 0 : else { return String(id + " must be a bool\n"); }
1469 : }
1470 119 : else{ return String(""); }
1471 : }
1472 :
1473 : // Read Vector<Int> from Record
1474 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
1475 : {
1476 0 : if( rec.isDefined( id ) )
1477 : {
1478 0 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1479 0 : { rec.get( id , val ); return String(""); }
1480 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1481 : {
1482 0 : Vector<Bool> vec; rec.get( id, vec );
1483 0 : if( vec.nelements()==0 ){val.resize(0); return String("");}
1484 0 : else{ return String(id + " must be a vector of strings.\n"); }
1485 0 : }
1486 0 : else { return String(id + " must be a vector of integers\n"); }
1487 : }
1488 0 : else{ return String(""); }
1489 : }
1490 :
1491 : // Read Vector<Float> from Record
1492 51 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1493 : {
1494 51 : if( rec.isDefined( id ) )
1495 : {
1496 51 : if( rec.dataType( id )==TpArrayFloat )
1497 : {
1498 17 : 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 34 : else if ( rec.dataType( id ) ==TpArrayDouble )
1512 : {
1513 0 : Vector<Double> vec; rec.get( id, vec );
1514 0 : val.resize(vec.nelements());
1515 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1516 0 : return String("");
1517 0 : }
1518 34 : else if ( rec.dataType( id ) ==TpArrayInt )
1519 : {
1520 0 : Vector<Int> vec; rec.get( id, vec );
1521 0 : val.resize(vec.nelements());
1522 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1523 0 : return String("");
1524 0 : }
1525 34 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1526 : {
1527 34 : Vector<Bool> vec; rec.get( id, vec );
1528 34 : 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 34 : }
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 17 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1539 : {
1540 17 : if( rec.isDefined( id ) )
1541 : {
1542 17 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1543 17 : { 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 102 : 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 102 : if ( casacore::Quantity::read( qa, instr ) ) { return String(""); }
1564 0 : else { return String("Error in converting " + instr + " to a Quantity \n"); }
1565 : }
1566 :
1567 : // Convert String to MDirection
1568 : // UR : TODO : If J2000 not specified, make it still work.
1569 0 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1570 : {
1571 : try
1572 : {
1573 0 : String tmpRF, tmpRA, tmpDEC;
1574 0 : std::istringstream iss(instr);
1575 0 : iss >> tmpRF >> tmpRA >> tmpDEC;
1576 0 : if( tmpDEC.length() == 0 )// J2000 may not be specified.
1577 : {
1578 0 : tmpDEC = tmpRA;
1579 0 : tmpRA = tmpRF;
1580 0 : tmpRF = String("J2000");
1581 : }
1582 0 : casacore::Quantity tmpQRA;
1583 0 : casacore::Quantity tmpQDEC;
1584 0 : casacore::Quantity::read(tmpQRA, tmpRA);
1585 0 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1586 :
1587 0 : if(tmpQDEC.getFullUnit()==Unit("deg") && tmpDEC.contains(":")){
1588 0 : LogIO os( LogOrigin("SynthesisParams","stringToMDirection",WHERE) );
1589 : os << LogIO::WARN
1590 : << "You provided the Declination/Latitude value \""<< tmpDEC
1591 : << "\" which is understood to be in units of hours.\n"
1592 : << "If you meant degrees, please replace \":\" by \".\"."
1593 0 : << LogIO::POST;
1594 0 : }
1595 :
1596 : MDirection::Types theRF;
1597 0 : Bool status = MDirection::getType(theRF, tmpRF);
1598 0 : if (!status) {
1599 0 : throw AipsError();
1600 : }
1601 0 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1602 0 : return String("");
1603 0 : }
1604 0 : catch(AipsError &x)
1605 : {
1606 0 : return String("Error in converting '" + instr + "' to MDirection. Need format : 'J2000 19:59:28.500 +40.44.01.50'\n");
1607 0 : }
1608 : }
1609 :
1610 : // Read Quantity from a Record string
1611 85 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1612 : {
1613 85 : if( rec.isDefined( id ) )
1614 : {
1615 68 : if( rec.dataType( id )==TpString )
1616 68 : { 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 17 : else{ return String(""); }
1620 : }
1621 :
1622 : // Read MDirection from a Record string
1623 17 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1624 : {
1625 17 : if( rec.isDefined( id ) )
1626 : {
1627 0 : if( rec.dataType( id )==TpString )
1628 0 : { String valstr; rec.get( id , valstr ); return stringToMDirection(valstr, val); }
1629 0 : else { return String(id + " must be a string in the format : 'J2000 19:59:28.500 +40.44.01.50'\n"); }
1630 : }
1631 17 : else{ return String(""); }
1632 : }
1633 :
1634 : // Convert MDirection to String
1635 0 : String SynthesisParams::MDirectionToString(MDirection val) const
1636 : {
1637 0 : MVDirection mvpc( val.getAngle() );
1638 0 : MVAngle ra = mvpc.get()(0);
1639 0 : MVAngle dec = mvpc.get()(1);
1640 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1641 0 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1642 0 : }
1643 :
1644 : // Convert Quantity to String
1645 0 : String SynthesisParams::QuantityToString(Quantity val) const
1646 : {
1647 0 : std::ostringstream ss;
1648 : //use digits10 to ensure the conersions involve use full decimal digits guranteed to be
1649 : //correct plus extra digits to deal with least significant digits (or replace with
1650 : // max_digits10 when it is available)
1651 0 : ss.precision(std::numeric_limits<double>::digits10+2);
1652 0 : ss << val;
1653 0 : return ss.str();
1654 : // NOTE - 2017.10.04: It was found (CAS-10773) that we cannot use to_string for this as
1655 : // the decimal place is fixed to 6 digits.
1656 : //TT: change to C++11 to_string which handles double value to string conversion
1657 : //return String(std::to_string( val.getValue(val.getUnit()) )) + val.getUnit() ;
1658 0 : }
1659 :
1660 : // Convert Record contains Quantity or Measure quantities to String
1661 0 : String SynthesisParams::recordQMToString(const Record &rec) const
1662 : {
1663 : Double val;
1664 0 : String unit;
1665 0 : if ( rec.isDefined("m0") )
1666 : {
1667 0 : Record subrec = rec.subRecord("m0");
1668 0 : subrec.get("value",val);
1669 0 : subrec.get("unit",unit);
1670 0 : }
1671 0 : else if (rec.isDefined("value") )
1672 : {
1673 0 : rec.get("value",val);
1674 0 : rec.get("unit",unit);
1675 : }
1676 0 : return String::toString(val) + unit;
1677 0 : }
1678 :
1679 :
1680 : /////////////////////// Selection Parameters
1681 :
1682 51 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1683 : {
1684 51 : setDefaults();
1685 51 : }
1686 :
1687 51 : SynthesisParamsSelect::~SynthesisParamsSelect()
1688 : {
1689 51 : }
1690 :
1691 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1692 0 : operator=(other);
1693 0 : }
1694 :
1695 17 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1696 17 : if(this!=&other) {
1697 17 : msname=other.msname;
1698 17 : spw=other.spw;
1699 17 : freqbeg=other.freqbeg;
1700 17 : freqend=other.freqend;
1701 17 : freqframe=other.freqframe;
1702 17 : field=other.field;
1703 17 : antenna=other.antenna;
1704 17 : timestr=other.timestr;
1705 17 : scan=other.scan;
1706 17 : obs=other.obs;
1707 17 : state=other.state;
1708 17 : uvdist=other.uvdist;
1709 17 : taql=other.taql;
1710 17 : usescratch=other.usescratch;
1711 17 : readonly=other.readonly;
1712 17 : incrmodel=other.incrmodel;
1713 17 : datacolumn=other.datacolumn;
1714 :
1715 : }
1716 17 : return *this;
1717 : }
1718 :
1719 34 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1720 : {
1721 34 : setDefaults();
1722 :
1723 34 : String err("");
1724 :
1725 : try
1726 : {
1727 :
1728 34 : err += readVal( inrec, String("msname"), msname );
1729 :
1730 34 : err += readVal( inrec, String("readonly"), readonly );
1731 34 : err += readVal( inrec, String("usescratch"), usescratch );
1732 :
1733 : // override with entries from savemodel.
1734 34 : String savemodel("");
1735 34 : err += readVal( inrec, String("savemodel"), savemodel );
1736 34 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1737 17 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1738 17 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1739 :
1740 34 : err += readVal( inrec, String("incrmodel"), incrmodel );
1741 :
1742 34 : err += readVal( inrec, String("spw"), spw );
1743 :
1744 : /// -------------------------------------------------------------------------------------
1745 : // Why are these params here ? Repeats of defineimage.
1746 34 : err += readVal( inrec, String("freqbeg"), freqbeg);
1747 34 : err += readVal( inrec, String("freqend"), freqend);
1748 :
1749 34 : String freqframestr( MFrequency::showType(freqframe) );
1750 34 : err += readVal( inrec, String("outframe"), freqframestr);
1751 34 : if( ! MFrequency::getType(freqframe, freqframestr) )
1752 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1753 : /// -------------------------------------------------------------------------------------
1754 :
1755 34 : err += readVal( inrec, String("field"),field);
1756 34 : err += readVal( inrec, String("antenna"),antenna);
1757 34 : err += readVal( inrec, String("timestr"),timestr);
1758 34 : err += readVal( inrec, String("scan"),scan);
1759 34 : err += readVal( inrec, String("obs"),obs);
1760 34 : err += readVal( inrec, String("state"),state);
1761 34 : err += readVal( inrec, String("uvdist"),uvdist);
1762 34 : err += readVal( inrec, String("taql"),taql);
1763 :
1764 34 : err += readVal( inrec, String("datacolumn"),datacolumn);
1765 :
1766 34 : err += verify();
1767 :
1768 34 : }
1769 0 : catch(AipsError &x)
1770 : {
1771 0 : err = err + x.getMesg() + "\n";
1772 0 : }
1773 :
1774 34 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1775 :
1776 34 : }
1777 :
1778 34 : String SynthesisParamsSelect::verify() const
1779 : {
1780 34 : String err;
1781 : // Does the MS exist on disk.
1782 34 : Directory thems( msname );
1783 34 : if( thems.exists() )
1784 : {
1785 : // Is it readable ?
1786 34 : if( ! thems.isReadable() )
1787 0 : { err += "MS " + msname + " is not readable.\n"; }
1788 : // Depending on 'readonly', is the MS writable ?
1789 34 : 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 68 : return err;
1796 34 : }
1797 :
1798 85 : void SynthesisParamsSelect::setDefaults()
1799 : {
1800 85 : msname="";
1801 85 : spw="";
1802 85 : freqbeg="";
1803 85 : freqend="";
1804 85 : MFrequency::getType(freqframe,"LSRK");
1805 85 : field="";
1806 85 : antenna="";
1807 85 : timestr="";
1808 85 : scan="";
1809 85 : obs="";
1810 85 : state="";
1811 85 : uvdist="";
1812 85 : taql="";
1813 85 : usescratch=false;
1814 85 : readonly=true;
1815 85 : incrmodel=false;
1816 85 : datacolumn="corrected";
1817 85 : }
1818 :
1819 17 : Record SynthesisParamsSelect::toRecord()const
1820 : {
1821 17 : Record selpar;
1822 17 : selpar.define("msname",msname);
1823 17 : selpar.define("spw",spw);
1824 17 : selpar.define("freqbeg",freqbeg);
1825 17 : selpar.define("freqend",freqend);
1826 17 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1827 : //looks like fromRecord looks for outframe !
1828 17 : selpar.define("outframe", MFrequency::showType(freqframe));
1829 17 : selpar.define("field",field);
1830 17 : selpar.define("antenna",antenna);
1831 17 : selpar.define("timestr",timestr);
1832 17 : selpar.define("scan",scan);
1833 17 : selpar.define("obs",obs);
1834 17 : selpar.define("state",state);
1835 17 : selpar.define("uvdist",uvdist);
1836 17 : selpar.define("taql",taql);
1837 17 : selpar.define("usescratch",usescratch);
1838 17 : selpar.define("readonly",readonly);
1839 17 : selpar.define("incrmodel",incrmodel);
1840 17 : selpar.define("datacolumn",datacolumn);
1841 :
1842 17 : return selpar;
1843 0 : }
1844 :
1845 :
1846 : /////////////////////// Image Parameters
1847 :
1848 51 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1849 : {
1850 51 : setDefaults();
1851 51 : }
1852 :
1853 51 : SynthesisParamsImage::~SynthesisParamsImage()
1854 : {
1855 51 : }
1856 :
1857 34 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1858 34 : if(this != &other){
1859 34 : imageName=other.imageName;
1860 34 : stokes=other.stokes;
1861 34 : startModel.resize(); startModel=other.startModel;
1862 34 : imsize.resize(); imsize=other.imsize;
1863 34 : cellsize.resize(); cellsize=other.cellsize;
1864 34 : projection=other.projection;
1865 34 : useNCP=other.useNCP;
1866 34 : phaseCenter=other.phaseCenter;
1867 34 : phaseCenterFieldId=other.phaseCenterFieldId;
1868 34 : obslocation=other.obslocation;
1869 34 : pseudoi=other.pseudoi;
1870 34 : nchan=other.nchan;
1871 34 : nTaylorTerms=other.nTaylorTerms;
1872 34 : chanStart=other.chanStart;
1873 34 : chanStep=other.chanStep;
1874 34 : freqStart=other.freqStart;
1875 34 : freqStep=other.freqStep;
1876 34 : refFreq=other.refFreq;
1877 34 : velStart=other.velStart;
1878 34 : velStep=other.velStep;
1879 34 : freqFrame=other.freqFrame;
1880 34 : mFreqStart=other.mFreqStart;
1881 34 : mFreqStep=other.mFreqStep;
1882 34 : mVelStart=other.mVelStart;
1883 34 : mVelStep=other.mVelStep;
1884 34 : restFreq.resize(); restFreq=other.restFreq;
1885 34 : start=other.start;
1886 34 : step=other.step;
1887 34 : frame=other.frame;
1888 34 : veltype=other.veltype;
1889 34 : mode=other.mode;
1890 34 : reffreq=other.reffreq;
1891 34 : sysvel=other.sysvel;
1892 34 : sysvelframe=other.sysvelframe;
1893 34 : sysvelvalue=other.sysvelvalue;
1894 34 : qmframe=other.qmframe;
1895 34 : mveltype=other.mveltype;
1896 34 : tststr=other.tststr;
1897 34 : startRecord=other.startRecord;
1898 34 : stepRecord=other.stepRecord;
1899 34 : reffreqRecord=other.reffreqRecord;
1900 34 : sysvelRecord=other.sysvelRecord;
1901 34 : restfreqRecord=other.restfreqRecord;
1902 34 : csysRecord=other.csysRecord;
1903 34 : csys=other.csys;
1904 34 : imshape.resize(); imshape=other.imshape;
1905 34 : freqFrameValid=other.freqFrameValid;
1906 34 : overwrite=other.overwrite;
1907 34 : deconvolver=other.deconvolver;
1908 34 : distance=other.distance;
1909 34 : trackDir=other.trackDir;
1910 34 : trackSource=other.trackSource;
1911 34 : movingSource=other.movingSource;
1912 :
1913 :
1914 :
1915 : }
1916 :
1917 34 : return *this;
1918 :
1919 : }
1920 :
1921 17 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec) {
1922 17 : fromRecord(inrec, False);
1923 17 : }
1924 :
1925 17 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec,
1926 : const casacore::Bool isSingleDish) {
1927 17 : setDefaults();
1928 17 : String err("");
1929 34 : LogIO os(LogOrigin("SynthesisParamsImage", "fromRecord"));
1930 :
1931 : try {
1932 17 : err += readVal(inrec, String("imagename"), imageName);
1933 :
1934 : //// imsize
1935 17 : if (inrec.isDefined("imsize")) {
1936 17 : DataType tp = inrec.dataType("imsize");
1937 17 : if (tp == TpInt) { // A single integer for both dimensions.
1938 : Int npix;
1939 0 : inrec.get("imsize", npix);
1940 0 : imsize.resize(2);
1941 0 : imsize.set(npix);
1942 17 : } else if (tp == TpArrayInt) { // An integer array : [ nx ] or [ nx, ny ]
1943 17 : Vector<Int> ims;
1944 17 : inrec.get("imsize", ims);
1945 17 : if (ims.nelements() == 1) { // [ nx ]
1946 0 : imsize.set(ims[0]);
1947 17 : } else if (ims.nelements() == 2) { // [ nx, ny ]
1948 17 : imsize[0] = ims[0];
1949 17 : 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 17 : } 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 17 : if (inrec.isDefined("cell")) {
1960 : try {
1961 17 : DataType tp = inrec.dataType("cell");
1962 17 : if (tp == TpInt ||
1963 17 : tp == TpFloat ||
1964 : tp == TpDouble) {
1965 0 : Double cell = inrec.asDouble("cell");
1966 0 : cellsize.set(Quantity(cell, "arcsec"));
1967 17 : } else if (tp == TpArrayInt ||
1968 17 : 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 17 : } else if (tp == TpString) {
1981 0 : String cell;
1982 0 : inrec.get("cell", cell);
1983 0 : Quantity qcell;
1984 0 : err += stringToQuantity(cell, qcell);
1985 0 : cellsize.set(qcell);
1986 17 : } else if (tp == TpArrayString) {
1987 17 : Array<String> cells;
1988 17 : inrec.get("cell", cells);
1989 17 : Vector<String> vcells(cells);
1990 17 : if (cells.nelements() == 1) { // [ cellx ]
1991 0 : Quantity qcell;
1992 0 : err += stringToQuantity(vcells[0], qcell);
1993 0 : cellsize.set(qcell);
1994 17 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1995 17 : err += stringToQuantity(vcells[0], cellsize[0]);
1996 17 : 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 17 : } 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 17 : err += readVal(inrec, String("stokes"), stokes);
2010 17 : if (stokes.matches("pseudoI")) {
2011 0 : stokes = "I";
2012 0 : pseudoi = true;
2013 : } else {
2014 17 : pseudoi = false;
2015 : }
2016 :
2017 : /// PseudoI
2018 :
2019 : ////nchan
2020 17 : err += readVal(inrec, String("nchan"), nchan);
2021 :
2022 : /// phaseCenter (as a string) . // Add INT support later.
2023 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2024 17 : if (inrec.isDefined("phasecenter")) {
2025 17 : String pcent("");
2026 17 : if (inrec.dataType("phasecenter") == TpString) {
2027 17 : inrec.get("phasecenter", pcent);
2028 17 : if (pcent.length() > 0) { // if it's zero length, it means 'figure out from first field in MS'.
2029 0 : err += readVal(inrec, String("phasecenter"), phaseCenter);
2030 0 : phaseCenterFieldId = -1;
2031 : //// Phase Center Field ID.... if explicitly specified, and not via phasecenter.
2032 : // Need this, to deal with a null phase center being translated to a string to go back out.
2033 0 : err += readVal(inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2034 : } else {
2035 17 : phaseCenterFieldId = -2; // deal with this later in buildCoordinateSystem to assign the first selected field
2036 : }
2037 : }
2038 17 : if (inrec.dataType("phasecenter") == TpInt ||
2039 51 : inrec.dataType("phasecenter") == TpFloat ||
2040 34 : 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 17 : if ((inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter") != TpInt
2045 17 : && inrec.dataType("phasecenter") != TpFloat && inrec.dataType("phasecenter") != TpDouble)) {
2046 0 : err += String("Cannot set phasecenter. Please specify a string or int\n");
2047 : }
2048 17 : }
2049 :
2050 : //// Projection
2051 17 : if (inrec.isDefined("projection")) {
2052 17 : if (inrec.dataType("projection") == TpString) {
2053 17 : String pstr;
2054 17 : inrec.get("projection", pstr);
2055 : try {
2056 17 : if (pstr.matches("NCP")) {
2057 0 : pstr = "SIN";
2058 0 : useNCP = true;
2059 : }
2060 17 : projection = Projection::type(pstr);
2061 0 : } catch (AipsError &x) {
2062 0 : err += String("Invalid projection code : " + pstr);
2063 0 : }
2064 17 : } else {
2065 0 : err += "projection must be a string\n";
2066 : }
2067 : } //projection
2068 :
2069 : // Frequency frame stuff.
2070 17 : err += readVal(inrec, String("specmode"), mode);
2071 : // Alias for 'mfs' is 'cont'
2072 17 : if (mode == "cont") mode = "mfs";
2073 :
2074 17 : err += readVal(inrec, String("outframe"), frame);
2075 17 : qmframe = "";
2076 : // mveltype is only set when start/step is given in mdoppler
2077 17 : mveltype = "";
2078 : //start
2079 17 : String startType("");
2080 17 : String widthType("");
2081 17 : if (inrec.isDefined("start")) {
2082 17 : if (inrec.dataType("start") == TpInt) {
2083 0 : err += readVal(inrec, String("start"), chanStart);
2084 0 : start = String::toString(chanStart);
2085 0 : startType = "chan";
2086 17 : } else if (inrec.dataType("start") == TpString) {
2087 17 : err += readVal(inrec, String("start"), start);
2088 17 : if (start.contains("Hz")) {
2089 0 : stringToQuantity(start, freqStart);
2090 0 : startType = "freq";
2091 17 : } 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 17 : if (inrec.isDefined("width")) {
2152 17 : if (inrec.dataType("width") == TpInt) {
2153 0 : err += readVal(inrec, String("width"), chanStep);
2154 0 : step = String::toString(chanStep);
2155 0 : widthType = "chan";
2156 17 : } else if (inrec.dataType("width") == TpString) {
2157 17 : err += readVal(inrec, String("width"), step);
2158 17 : if (step.contains("Hz")) {
2159 0 : stringToQuantity(step, freqStep);
2160 0 : widthType = "freq";
2161 17 : } 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 17 : 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 17 : if (inrec.isDefined("reffreq")) {
2226 17 : if (inrec.dataType("reffreq") == TpString) {
2227 17 : 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 17 : err += readVal(inrec, String("veltype"), veltype);
2254 17 : veltype = mveltype != "" ? mveltype : veltype;
2255 : // sysvel (String, Quantity)
2256 17 : if (inrec.isDefined("sysvel")) {
2257 17 : if (inrec.dataType("sysvel") == TpString) {
2258 17 : 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 17 : err += readVal(inrec, String("sysvelframe"), sysvelframe);
2270 :
2271 : // rest frequencies (record or vector of Strings)
2272 17 : if (inrec.isDefined("restfreq")) {
2273 17 : Vector<String> rfreqs(0);
2274 17 : Record restfreqSubRecord;
2275 17 : 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 17 : } else if (inrec.dataType("restfreq") == TpArrayString) {
2286 0 : err += readVal(inrec, String("restfreq"), rfreqs);
2287 17 : } else if (inrec.dataType("restfreq") == TpString) {
2288 0 : rfreqs.resize(1);
2289 0 : err += readVal(inrec, String("restfreq"), rfreqs[0]);
2290 : }
2291 17 : restFreq.resize(rfreqs.nelements());
2292 17 : for (uInt fr = 0; fr < rfreqs.nelements(); fr++) {
2293 0 : err += stringToQuantity(rfreqs[fr], restFreq[fr]);
2294 : }
2295 17 : } // if def restfreq
2296 :
2297 : // optional - coordsys, imshape
2298 : // if exist use them. May need a consistency check with the rest of impars?
2299 17 : if (inrec.isDefined("csys")) {
2300 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
2301 0 : if (inrec.dataType("csys") == TpRecord) {
2302 : //csysRecord = inrec.subRecord("csys");
2303 0 : csysRecord.defineRecord("coordsys", inrec.subRecord("csys"));
2304 : }
2305 0 : if (inrec.isDefined("imshape")) {
2306 0 : if (inrec.dataType("imshape") == TpArrayInt) {
2307 0 : err += readVal(inrec, String("imshape"), imshape);
2308 : }
2309 : }
2310 : }
2311 :
2312 : //String freqframestr( MFrequency::showType(freqFrame) );
2313 : //err += readVal( inrec, String("outframe"), freqframestr);
2314 : //if( ! MFrequency::getType(freqFrame, freqframestr) )
2315 : // { err += "Invalid Frequency Frame " + freqframestr ; }
2316 :
2317 17 : String freqframestr = (frame != "" && qmframe != "") ? qmframe : frame;
2318 17 : if (frame != "" && !MFrequency::getType(freqFrame, freqframestr)) {
2319 0 : err += "Invalid Frequency Frame " + freqframestr;
2320 : }
2321 17 : err += readVal(inrec, String("restart"), overwrite);
2322 :
2323 17 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2324 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2325 17 : if (inrec.isDefined("startmodel")) {
2326 17 : if (inrec.dataType("startmodel") == TpString) {
2327 17 : String onemodel;
2328 17 : err += readVal(inrec, String("startmodel"), onemodel);
2329 17 : if (onemodel.length() > 0) {
2330 0 : startModel.resize(1);
2331 0 : startModel[0] = onemodel;
2332 : } else {
2333 17 : startModel.resize();
2334 : }
2335 17 : } else if (inrec.dataType("startmodel") == TpArrayString ||
2336 0 : inrec.dataType("startmodel") == TpArrayBool) {
2337 0 : err += readVal(inrec, String("startmodel"), startModel);
2338 : } else {
2339 0 : err += String("startmodel must be either a string(singleterm) or a list of strings(multiterm)\n");
2340 : }
2341 : }
2342 :
2343 17 : err += readVal(inrec, String("nterms"), nTaylorTerms);
2344 17 : err += readVal(inrec, String("deconvolver"), deconvolver);
2345 :
2346 : // Force nchan=1 for anything other than cube modes...
2347 17 : if (mode == "mfs") nchan = 1;
2348 : //read obslocation
2349 17 : if (inrec.isDefined("obslocation_rec")) {
2350 0 : String errorobs;
2351 0 : const Record obsrec = inrec.asRecord("obslocation_rec");
2352 0 : MeasureHolder mh;
2353 0 : if (!mh.fromRecord(errorobs, obsrec)) {
2354 0 : err += errorobs;
2355 : }
2356 0 : obslocation = mh.asMPosition();
2357 0 : }
2358 17 : err += verify(isSingleDish);
2359 17 : }
2360 0 : catch (AipsError &x) {
2361 0 : err = err + x.getMesg() + "\n";
2362 0 : }
2363 :
2364 17 : if (err.length() > 0) throw(AipsError("Invalid Image Parameter set : " + err));
2365 17 : }
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 17 : String SynthesisParamsImage::verify(const casacore::Bool isSingleDish) const
2401 : {
2402 17 : LogIO os;
2403 17 : String err;
2404 :
2405 17 : if(imageName=="") {err += "Please supply an image name\n";}
2406 :
2407 17 : if(imsize.nelements() != 2){err += "imsize must be a vector of 2 Ints\n";}
2408 17 : if(cellsize.nelements() != 2) {err += "cellsize must be a vector of 2 Quantities\n";}
2409 17 : 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 17 : if((mode=="mfs") && nchan > 1){
2418 0 : err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";
2419 : }
2420 :
2421 17 : if(! stokes.matches("I") && ! stokes.matches("Q") &&
2422 0 : ! stokes.matches("U") && ! stokes.matches("V") &&
2423 0 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2424 0 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2425 0 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2426 0 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2427 0 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2428 17 : ! 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 17 : 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 17 : Int imxnew = SynthesisUtilMethods::getOptimumSize(imsize[0]);
2468 17 : Int imynew = SynthesisUtilMethods::getOptimumSize(imsize[1]);
2469 17 : 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 34 : return err;
2483 17 : }// 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 68 : void SynthesisParamsImage::setDefaults()
2496 : {
2497 : // Image definition parameters
2498 68 : imageName = String("");
2499 68 : imsize.resize(2); imsize.set(100);
2500 68 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2501 68 : stokes="I";
2502 68 : phaseCenter=MDirection();
2503 68 : phaseCenterFieldId=-1;
2504 68 : projection=Projection::SIN;
2505 68 : useNCP=false;
2506 68 : startModel=Vector<String>(0);
2507 68 : freqFrameValid=True;
2508 68 : overwrite=false;
2509 : // PseudoI
2510 68 : pseudoi=false;
2511 :
2512 : // Spectral coordinates
2513 68 : nchan=1;
2514 68 : mode="mfs";
2515 68 : start="";
2516 68 : step="";
2517 68 : chanStart=0;
2518 68 : 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 68 : freqStart=Quantity(0,"");
2524 68 : freqStep=Quantity(0,"");
2525 68 : velStart=Quantity(0,"");
2526 68 : velStep=Quantity(0,"");
2527 68 : veltype=String("radio");
2528 68 : restFreq.resize(0);
2529 68 : refFreq = Quantity(0,"Hz");
2530 68 : frame = "";
2531 68 : freqFrame=MFrequency::LSRK;
2532 68 : sysvel="";
2533 68 : sysvelframe="";
2534 68 : sysvelvalue=Quantity(0.0,"m/s");
2535 68 : nTaylorTerms=1;
2536 68 : deconvolver="hogbom";
2537 : ///csysRecord=Record();
2538 68 : }
2539 :
2540 0 : Record SynthesisParamsImage::toRecord() const
2541 : {
2542 0 : Record impar;
2543 0 : impar.define("imagename", imageName);
2544 0 : impar.define("imsize", imsize);
2545 0 : Vector<String> cells(2);
2546 0 : cells[0] = QuantityToString( cellsize[0] );
2547 0 : cells[1] = QuantityToString( cellsize[1] );
2548 0 : impar.define("cell", cells );
2549 0 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2550 0 : else{impar.define("stokes", stokes);}
2551 0 : impar.define("nchan", nchan);
2552 0 : impar.define("nterms", nTaylorTerms);
2553 0 : impar.define("deconvolver",deconvolver);
2554 0 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2555 0 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2556 0 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2557 :
2558 0 : impar.define("specmode", mode );
2559 : // start and step can be one of these types
2560 0 : if( start!="" )
2561 : {
2562 0 : if( !start.contains("Hz") && !start.contains("m/s") &&
2563 0 : String::toInt(start) == chanStart )
2564 : {
2565 0 : impar.define("start",chanStart);
2566 : }
2567 0 : else if( startRecord.nfields() > 0 )
2568 : {
2569 0 : impar.defineRecord("start", startRecord );
2570 : }
2571 : else
2572 : {
2573 0 : impar.define("start",start);
2574 : }
2575 : }
2576 : else {
2577 0 : impar.define("start", start );
2578 : }
2579 0 : if( step!="" )
2580 : {
2581 0 : if( !step.contains("Hz") && !step.contains("m/s") &&
2582 0 : String::toInt(step) == chanStep )
2583 : {
2584 0 : impar.define("width", chanStep);
2585 : }
2586 0 : else if( stepRecord.nfields() > 0 )
2587 : {
2588 0 : impar.defineRecord("width",stepRecord);
2589 : }
2590 : else
2591 : {
2592 0 : impar.define("width",step);
2593 : }
2594 : }
2595 : else
2596 : {
2597 0 : impar.define("width", step);
2598 : }
2599 : //impar.define("chanstart", chanStart );
2600 : //impar.define("chanstep", chanStep );
2601 : //impar.define("freqstart", QuantityToString( freqStart ));
2602 : //impar.define("freqstep", QuantityToString( freqStep ) );
2603 : //impar.define("velstart", QuantityToString( velStart ));
2604 : //impar.define("velstep", QuantityToString( velStep ) );
2605 0 : impar.define("veltype", veltype);
2606 0 : if (restfreqRecord.nfields() != 0 )
2607 : {
2608 0 : impar.defineRecord("restfreq", restfreqRecord);
2609 : }
2610 : else
2611 : {
2612 0 : Vector<String> rfs( restFreq.nelements() );
2613 0 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2614 0 : impar.define("restfreq", rfs);
2615 0 : }
2616 : //impar.define("reffreq", QuantityToString(refFreq));
2617 : //reffreq
2618 0 : if( reffreqRecord.nfields() != 0 )
2619 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2620 : else
2621 0 : { impar.define("reffreq",reffreq); }
2622 : //impar.define("reffreq", reffreq );
2623 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2624 0 : impar.define("outframe", frame );
2625 : //sysvel
2626 0 : if( sysvelRecord.nfields() != 0 )
2627 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2628 : else
2629 0 : { impar.define("sysvel", sysvel );}
2630 0 : impar.define("sysvelframe", sysvelframe );
2631 :
2632 0 : impar.define("restart",overwrite );
2633 0 : impar.define("freqframevalid", freqFrameValid);
2634 0 : impar.define("startmodel", startModel );
2635 :
2636 0 : if( csysRecord.isDefined("coordsys") )
2637 : {
2638 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2639 0 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2640 0 : impar.define("imshape", imshape);
2641 : }
2642 : // else cout << " NO CSYS INFO to write to output record " << endl;
2643 : ///Now save obslocation
2644 0 : Record tmprec;
2645 0 : String err;
2646 0 : MeasureHolder mh(obslocation);
2647 0 : if(mh.toRecord(err, tmprec)){
2648 0 : impar.defineRecord("obslocation_rec", tmprec);
2649 : }
2650 : else{
2651 0 : throw(AipsError("failed to save obslocation to record"));
2652 : }
2653 0 : return impar;
2654 0 : }
2655 :
2656 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2657 : /////////////////////////// Build a coordinate system. ////////////////////////////////////////
2658 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2659 : //// To be used from SynthesisImager, to construct the images it needs
2660 : //// To also be connected to a 'makeimage' method of the synthesisimager tool.
2661 : //// ( need to supply MS only to add 'ObsInfo' to the csys )
2662 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2663 :
2664 :
2665 :
2666 17 : 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 17 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2674 17 : vi2.originChunks();
2675 17 : vi2.origin();
2676 : /// This version uses the new vi2/vb2
2677 : // get the first ms for multiple MSes
2678 : //MeasurementSet msobj=vi2.ms();
2679 17 : Int fld=vb->fieldId()(0);
2680 :
2681 : //handling first ms only
2682 17 : Double gfreqmax=-1.0;
2683 17 : Double gdatafend=-1.0;
2684 17 : Double gdatafstart=1e14;
2685 17 : Double gfreqmin=1e14;
2686 17 : Vector<Int> spwids0;
2687 17 : Int j=0;
2688 17 : Int minfmsid=0;
2689 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2690 17 : Double imStartFreq=getCubeImageStartFreq();
2691 17 : std::vector<Int> sourceMsWithStartFreq;
2692 :
2693 :
2694 34 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2695 : //auto forMS0=chansel.find(0);
2696 17 : map<Int, Vector<Int> > spwsels=forMS0->second;
2697 17 : Int nspws=spwsels.size();
2698 17 : Vector<Int> spwids(nspws);
2699 17 : Vector<Int> nChannels(nspws);
2700 17 : Vector<Int> firstChannels(nspws);
2701 : //Vector<Int> channelIncrement(nspws);
2702 :
2703 17 : Int k=0;
2704 34 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2705 17 : spwids[k]=it->first;
2706 17 : nChannels[k]=(it->second)[0];
2707 17 : firstChannels[k]=(it->second)[1];
2708 : }
2709 17 : if(j==0) {
2710 17 : spwids0.resize();
2711 17 : 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 17 : Double freqmin=0, freqmax=0;
2723 17 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2724 :
2725 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2726 17 : 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 17 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2732 : //cerr << "after " << datafstart << " " << datafend << endl;
2733 17 : if((datafstart > datafend) || !status)
2734 0 : throw(AipsError("spw selection failed"));
2735 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2736 :
2737 17 : if (mode=="cubedata") {
2738 0 : freqmin = datafstart;
2739 0 : freqmax = datafend;
2740 : }
2741 17 : 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 34 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2766 34 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2767 : //cerr << "after " << freqmin << " " << freqmax << endl;
2768 : }
2769 :
2770 :
2771 :
2772 :
2773 17 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2774 17 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2775 17 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2776 17 : 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 17 : 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 17 : }
2802 17 : 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 17 : MeasurementSet msobj = *mss[minfmsid];
2809 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2810 34 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2811 17 : }
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 17 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2844 : MeasurementSet& msobj,
2845 : Vector<Int> spwids, Int fld,
2846 : Double freqmin, Double freqmax,
2847 : Double datafstart, Double datafend )
2848 : {
2849 34 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2850 :
2851 17 : CoordinateSystem csys;
2852 17 : 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 17 : MSColumns msc(msobj);
2865 17 : String telescop = msc.observation().telescopeName()(0);
2866 17 : MEpoch obsEpoch = msc.timeMeas()(0);
2867 17 : MPosition obsPosition;
2868 17 : 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 17 : MDirection phaseCenterToUse = phaseCenter;
2880 :
2881 17 : if( phaseCenterFieldId != -1 )
2882 : {
2883 17 : MSFieldColumns msfield(msobj.field());
2884 17 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2885 : {
2886 17 : if(trackSource){
2887 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2888 : }
2889 : else{
2890 17 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2891 : }
2892 : }
2893 : else
2894 : {
2895 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2896 : }
2897 17 : }
2898 : // Setup Phase center if it is specified only by field id.
2899 :
2900 : /////////////////// Direction Coordinates
2901 17 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2902 : // Normalize correctly
2903 17 : MVAngle ra=mvPhaseCenter.get()(0);
2904 17 : ra(0.0);
2905 :
2906 17 : MVAngle dec=mvPhaseCenter.get()(1);
2907 17 : Vector<Double> refCoord(2);
2908 17 : refCoord(0)=ra.get().getValue();
2909 17 : refCoord(1)=dec;
2910 17 : Vector<Double> refPixel(2);
2911 17 : refPixel(0) = Double(imsize[0]/2);
2912 17 : refPixel(1) = Double(imsize[1]/2);
2913 :
2914 17 : Vector<Double> deltas(2);
2915 17 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2916 17 : deltas(1)=cellsize[1].get("rad").getValue();
2917 17 : Matrix<Double> xform(2,2);
2918 17 : xform=0.0;xform.diagonal()=1.0;
2919 :
2920 17 : Vector<Double> projparams(2);
2921 17 : projparams = 0.0;
2922 17 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2923 17 : Projection projTo( projection.type(), projparams );
2924 :
2925 : DirectionCoordinate
2926 17 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2927 : // projection,
2928 : projTo,
2929 51 : refCoord(0), refCoord(1),
2930 51 : deltas(0), deltas(1),
2931 : xform,
2932 34 : 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 17 : obslocation=obsPosition;
2940 17 : ObsInfo myobsinfo;
2941 17 : myobsinfo.setTelescope(telescop);
2942 17 : myobsinfo.setPointingCenter(mvPhaseCenter);
2943 17 : myobsinfo.setObsDate(obsEpoch);
2944 17 : 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 17 : if(spwids.nelements()==0)
2960 : {
2961 0 : Int nspw=msc.spectralWindow().nrow();
2962 0 : spwids.resize(nspw);
2963 0 : indgen(spwids);
2964 : }
2965 17 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2966 17 : Vector<Double> dataChanFreq, dataChanWidth;
2967 17 : std::vector<std::vector<Int> > averageWhichChan;
2968 17 : std::vector<std::vector<Int> > averageWhichSPW;
2969 17 : std::vector<std::vector<Double> > averageChanFrac;
2970 :
2971 17 : if(spwids.nelements()==1)
2972 : {
2973 17 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2974 17 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
2975 : }
2976 : else
2977 : {
2978 0 : if(!MSTransformRegridder::combineSpwsCore(os,msobj, spwids,dataChanFreq,dataChanWidth,
2979 : averageWhichChan,averageWhichSPW,averageChanFrac))
2980 : {
2981 0 : os << LogIO::SEVERE << "Error combining SpWs" << LogIO::POST;
2982 : }
2983 : }
2984 17 : Double minDataFreq = min(dataChanFreq);
2985 17 : 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 17 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3037 17 : String cubemode;
3038 17 : if ( qrestfreq.getValue("Hz")==0 )
3039 : {
3040 17 : MSDopplerUtil msdoppler(msobj);
3041 17 : Vector<Double> restfreqvec;
3042 17 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3043 17 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3044 17 : 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 17 : }
3062 : Double refPix;
3063 17 : Vector<Double> chanFreq;
3064 17 : Vector<Double> chanFreqStep;
3065 17 : String specmode;
3066 :
3067 17 : 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 17 : 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 17 : Bool nonLinearFreq(false);
3083 17 : String veltype_p=veltype;
3084 17 : veltype_p.upcase();
3085 51 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3086 51 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3087 : {
3088 0 : nonLinearFreq= true;
3089 : }
3090 :
3091 17 : SpectralCoordinate mySpectral;
3092 : Double stepf;
3093 17 : if(!nonLinearFreq)
3094 : //TODO: velocity mode default start case (use last channels?)
3095 : {
3096 17 : Double startf=chanFreq[0];
3097 : //Double stepf=chanFreqStep[0];
3098 17 : if(chanFreq.nelements()==1)
3099 : {
3100 17 : stepf=chanFreqStep[0];
3101 : }
3102 : else
3103 : {
3104 0 : stepf=chanFreq[1]-chanFreq[0];
3105 : }
3106 17 : Double restf=qrestfreq.getValue("Hz");
3107 : //stepf=9e8;
3108 17 : 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 17 : 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 17 : 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 34 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3129 17 : 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 17 : uInt nrestfreq = restFreq.nelements();
3159 17 : 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 17 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3177 17 : if(whichStokes.nelements()==0)
3178 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3179 17 : StokesCoordinate myStokes(whichStokes);
3180 :
3181 : ////////////////// Build Full coordinate system.
3182 :
3183 : //CoordinateSystem csys;
3184 17 : csys.addCoordinate(myRaDec);
3185 17 : csys.addCoordinate(myStokes);
3186 17 : csys.addCoordinate(mySpectral);
3187 17 : csys.setObsInfo(myobsinfo);
3188 :
3189 : //store back csys to impars record
3190 : //cerr<<"save csys to csysRecord..."<<endl;
3191 17 : if(csysRecord.isDefined("coordsys"))
3192 0 : csysRecord.removeField("coordsys");
3193 17 : csys.save(csysRecord,"coordsys");
3194 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3195 : // imshape
3196 17 : imshape.resize(4);
3197 17 : imshape[0] = imsize[0];
3198 17 : imshape[1] = imsize[0];
3199 17 : imshape[2] = whichStokes.nelements();
3200 17 : 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 17 : } // end of else when coordsys record is not defined...
3234 :
3235 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3236 :
3237 34 : return csys;
3238 17 : }
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 17 : 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 17 : String inStart, inStep;
3531 17 : specmode = findSpecMode(mode);
3532 17 : String freqframe;
3533 17 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3534 34 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3535 :
3536 17 : refPix=0.0;
3537 17 : Bool descendingfreq(false);
3538 17 : Bool descendingoutfreq(false);
3539 :
3540 17 : if( mode.contains("cube") )
3541 : {
3542 0 : String restfreq=QuantityToString(qrestfreq);
3543 : // use frame from input start or width in MFreaquency or MRadialVelocity
3544 0 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3545 : // emit warning here if qmframe is used
3546 : //
3547 0 : inStart = start;
3548 0 : inStep = step;
3549 0 : if( specmode=="channel" )
3550 : {
3551 0 : inStart = String::toString(chanStart);
3552 0 : inStep = String::toString(chanStep);
3553 : // negative step -> descending channel indices
3554 0 : if (inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3555 : // input frame is the data frame
3556 : //freqframe = MFrequency::showType(dataFrame);
3557 : }
3558 0 : else if( specmode=="frequency" )
3559 : {
3560 : //if ( freqStart.getValue("Hz") == 0 && freqStart.getUnit() != "" ) { // default start
3561 : //start = String::toString( freqmin ) + freqStart.getUnit();
3562 : //}
3563 : //else {
3564 : //start = String::toString( freqStart.getValue(freqStart.getUnit()) )+freqStart.getUnit();
3565 : //}
3566 : //step = String::toString( freqStep.getValue(freqStep.getUnit()) )+freqStep.getUnit();
3567 : // negative freq width -> descending freq ordering
3568 0 : if(inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3569 : }
3570 0 : else if( specmode=="velocity" )
3571 : {
3572 : // if velStart is empty set start to vel of freqmin or freqmax?
3573 : //if ( velStart.getValue(velStart.getUnit()) == 0 && !(velStart.getUnit().contains("m/s")) ) {
3574 : // start = "";
3575 : //}
3576 : //else {
3577 : // start = String::toString( velStart.getValue(velStart.getUnit()) )+velStart.getUnit();
3578 : //}
3579 : //step = String::toString( velStep.getValue(velStep.getUnit()) )+velStep.getUnit();
3580 : // positive velocity width -> descending freq ordering
3581 0 : if (!inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3582 : }
3583 :
3584 0 : if (inStep=='0') inStep="";
3585 :
3586 0 : MRadialVelocity mSysVel;
3587 0 : Quantity qVel;
3588 : MRadialVelocity::Types mRef;
3589 0 : if(mode!="cubesource")
3590 : {
3591 :
3592 :
3593 0 : if(freqframe=="SOURCE")
3594 : {
3595 : os << LogIO::SEVERE << "freqframe=\"SOURCE\" is only allowed for mode=\"cubesrc\""
3596 0 : << LogIO::EXCEPTION;
3597 0 : return false;
3598 : }
3599 : }
3600 : else // only for cubesrc mode: TODO- check for the ephemeris info.
3601 : {
3602 0 : freqframe=MFrequency::showType(dataFrame);
3603 0 : if(sysvel!="") {
3604 0 : stringToQuantity(sysvel,qVel);
3605 0 : MRadialVelocity::getType(mRef,sysvelframe);
3606 0 : mSysVel=MRadialVelocity(qVel,mRef);
3607 : }
3608 : else // and if no ephemeris info, issue a warning...
3609 0 : { mSysVel=MRadialVelocity();}
3610 : }
3611 : // cubedata mode: input start, step are those of the input data frame
3612 0 : if ( mode=="cubedata" )
3613 : {
3614 0 : freqframe=MFrequency::showType(dataFrame);
3615 0 : freqFrameValid=false; // no conversion for vb.lsrfrequency()
3616 : }
3617 : //if ( mode=="cubedata" ) freqframe=MFrequency::REST;
3618 :
3619 : // *** NOTE ***
3620 : // calcChanFreqs alway returns chanFreq in
3621 : // ascending freq order.
3622 : // for step < 0 calcChanFreqs returns chanFreq that
3623 : // contains start freq. in its last element of the vector.
3624 : //
3625 0 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3626 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3627 0 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3628 0 : << LogIO::POST;
3629 0 : ostringstream ostr;
3630 0 : ostr << " phaseCenter='" << phaseCenter;
3631 0 : os << String(ostr)<<"' ";
3632 :
3633 : Double dummy; // dummy variable - weightScale is not used here
3634 0 : Bool rst=MSTransformRegridder::calcChanFreqs(os,
3635 : chanFreq,
3636 : chanFreqStep,
3637 : dummy,
3638 : dataChanFreq,
3639 : dataChanWidth,
3640 : phaseCenter,
3641 : dataFrame,
3642 : obsEpoch,
3643 : obsPosition,
3644 : specmode,
3645 : nchan,
3646 : inStart,
3647 : inStep,
3648 : restfreq,
3649 : freqframe,
3650 0 : veltype,
3651 : verbose,
3652 : mSysVel
3653 : );
3654 :
3655 0 : if( nchan==-1 )
3656 : {
3657 0 : nchan = chanFreq.nelements();
3658 0 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3659 : }
3660 :
3661 0 : if (!rst) {
3662 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
3663 0 : << LogIO::EXCEPTION;
3664 0 : return false;
3665 : }
3666 : os << LogIO::DEBUG1
3667 0 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3668 0 : << LogIO::POST;
3669 :
3670 0 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3671 0 : descendingoutfreq = true;
3672 : }
3673 :
3674 : //if (descendingfreq && !descendingoutfreq) {
3675 0 : if ((specmode=="channel" && descendingfreq==1)
3676 0 : || (specmode!="channel" && (descendingfreq != descendingoutfreq))) {
3677 : // reverse the freq vector if necessary so the first element can be
3678 : // used to set spectralCoordinates in all the cases.
3679 : //
3680 : // also do for chanFreqStep..
3681 0 : std::vector<Double> stlchanfreq;
3682 0 : chanFreq.tovector(stlchanfreq);
3683 0 : std::reverse(stlchanfreq.begin(),stlchanfreq.end());
3684 0 : chanFreq=Vector<Double>(stlchanfreq);
3685 0 : chanFreqStep=-chanFreqStep;
3686 0 : }
3687 0 : }
3688 17 : else if ( mode=="mfs" ) {
3689 17 : chanFreq.resize(1);
3690 17 : chanFreqStep.resize(1);
3691 : //chanFreqStep[0] = freqmax - freqmin;
3692 17 : Double freqmean = (freqmin + freqmax)/2;
3693 17 : if (refFreq.getValue("Hz")==0) {
3694 17 : chanFreq[0] = freqmean;
3695 17 : refPix = 0.0;
3696 17 : chanFreqStep[0] = freqmax - freqmin;
3697 : }
3698 : else {
3699 0 : chanFreq[0] = refFreq.getValue("Hz");
3700 : // Set the new reffreq to be the refPix (CAS-9518)
3701 0 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3702 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3703 0 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3704 : }
3705 :
3706 17 : if( nchan==-1 ) nchan=1;
3707 17 : 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 17 : return true;
3719 :
3720 17 : }//getImFreq
3721 : /////////////////////////
3722 17 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3723 17 : Double inStartFreq=-1.0;
3724 17 : String checkspecmode("");
3725 17 : if(mode.contains("cube")) {
3726 0 : checkspecmode = findSpecMode(mode);
3727 : }
3728 17 : if(checkspecmode!="") {
3729 0 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3730 0 : if(checkspecmode=="channel") {
3731 0 : inStartFreq=-1.0;
3732 : }
3733 : else {
3734 0 : if(checkspecmode=="frequency") {
3735 0 : inStartFreq = freqStart.get("Hz").getValue();
3736 : }
3737 0 : else if(checkspecmode=="velocity") {
3738 : MDoppler::Types DopType;
3739 0 : MDoppler::getType(DopType, veltype);
3740 0 : MDoppler mdop(velStart,DopType);
3741 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3742 0 : inStartFreq = MFrequency::fromDoppler(mdop, qrestfreq.getValue(Unit("Hz")), mfreqframe).getValue();
3743 0 : }
3744 : }
3745 : }
3746 :
3747 17 : return inStartFreq;
3748 :
3749 17 : }
3750 :
3751 17 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3752 : {
3753 17 : String specmode;
3754 17 : specmode="channel";
3755 17 : if ( mode.contains("cube") ) {
3756 : // if velstart or velstep is defined -> specmode='vel'
3757 : // else if freqstart or freqstep is defined -> specmode='freq'
3758 : // velocity: assume unset if velStart => 0.0 with no unit
3759 : // assume unset if velStep => 0.0 with/without unit
3760 0 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3761 0 : !( velStep.getValue()==0.0)) {
3762 0 : specmode="velocity";
3763 : }
3764 0 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3765 0 : !(freqStep.getValue()==0.0)) {
3766 0 : specmode="frequency";
3767 : }
3768 : }
3769 17 : return specmode;
3770 0 : }
3771 :
3772 :
3773 51 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3774 : {
3775 51 : Vector<Int> whichStokes(0);
3776 51 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3777 0 : stokes=="RR" ||stokes=="LL" ||
3778 51 : stokes=="XX" || stokes=="YY" ) {
3779 51 : whichStokes.resize(1);
3780 51 : whichStokes(0)=Stokes::type(stokes);
3781 : }
3782 0 : else if(stokes=="IV" || stokes=="IQ" ||
3783 0 : stokes=="RRLL" || stokes=="XXYY" ||
3784 0 : stokes=="QU" || stokes=="UV"){
3785 0 : whichStokes.resize(2);
3786 :
3787 0 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3788 0 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3789 0 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3790 0 : else if(stokes=="XXYY"){whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
3791 0 : else if(stokes=="QU"){whichStokes[0]=Stokes::Q; whichStokes[1]=Stokes::U; }
3792 0 : else if(stokes=="UV"){ whichStokes[0]=Stokes::U; whichStokes[1]=Stokes::V; }
3793 :
3794 : }
3795 :
3796 0 : else if(stokes=="IQU" || stokes=="IUV") {
3797 0 : whichStokes.resize(3);
3798 0 : if(stokes=="IUV")
3799 0 : {whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::U; whichStokes[2]=Stokes::V;}
3800 : else
3801 0 : {whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q; whichStokes[2]=Stokes::U;}
3802 : }
3803 0 : else if(stokes=="IQUV"){
3804 0 : whichStokes.resize(4);
3805 0 : whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;
3806 0 : whichStokes(2)=Stokes::U; whichStokes(3)=Stokes::V;
3807 : }
3808 :
3809 51 : return whichStokes;
3810 0 : }// decidenpolplanes
3811 :
3812 34 : IPosition SynthesisParamsImage::shp() const
3813 : {
3814 34 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3815 :
3816 34 : 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 68 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3822 : }
3823 :
3824 17 : Record SynthesisParamsImage::getcsys() const
3825 : {
3826 17 : 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 51 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3883 : {
3884 51 : setDefaults();
3885 51 : }
3886 :
3887 51 : SynthesisParamsGrid::~SynthesisParamsGrid()
3888 : {
3889 51 : }
3890 :
3891 :
3892 17 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3893 : {
3894 17 : setDefaults();
3895 :
3896 17 : String err("");
3897 :
3898 : try {
3899 17 : err += readVal( inrec, String("imagename"), imageName );
3900 :
3901 : // FTMachine parameters
3902 17 : err += readVal( inrec, String("gridder"), gridder );
3903 17 : err += readVal( inrec, String("padding"), padding );
3904 17 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3905 17 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3906 17 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3907 17 : err += readVal( inrec, String("convfunc"), convFunc );
3908 :
3909 17 : err += readVal( inrec, String("vptable"), vpTable );
3910 :
3911 :
3912 :
3913 : // convert 'gridder' to 'ftmachine' and 'mtype'
3914 17 : ftmachine = "gridft";
3915 17 : mType = "default";
3916 17 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
3917 17 : ftmachine = "gridft";
3918 : }
3919 0 : else if ( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) &&
3920 0 : (wprojplanes>1 || wprojplanes==-1) ) {
3921 0 : ftmachine = "wprojectft";
3922 : }
3923 : //facetting alone use gridft
3924 0 : else if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes==1))
3925 0 : {ftmachine=="gridft";}
3926 :
3927 0 : else if (gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" ) {
3928 0 : ftmachine = "mosaicft";
3929 : }
3930 :
3931 0 : else if (gridder=="imagemosaic") {
3932 0 : mType = "imagemosaic";
3933 0 : if (wprojplanes>1 || wprojplanes==-1) {
3934 0 : ftmachine = "wprojectft";
3935 : }
3936 : }
3937 :
3938 0 : else if (gridder=="awproject" || gridder=="awprojectft" || gridder=="awp") {
3939 0 : ftmachine = "awprojectft";
3940 : }
3941 :
3942 0 : else if (gridder=="singledish") {
3943 :
3944 0 : ftmachine = "sd";
3945 : }
3946 : else{
3947 0 : ftmachine=gridder;
3948 0 : ftmachine.downcase();
3949 :
3950 : }
3951 :
3952 17 : String deconvolver;
3953 17 : err += readVal( inrec, String("deconvolver"), deconvolver );
3954 17 : if (deconvolver=="mtmfs") {
3955 0 : mType = "multiterm"; // Takes precedence over imagemosaic
3956 : }
3957 :
3958 : // facets
3959 17 : err += readVal( inrec, String("facets"), facets );
3960 : // chanchunks
3961 17 : err += readVal( inrec, String("chanchunks"), chanchunks );
3962 :
3963 : // Spectral interpolation
3964 17 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
3965 : // Track moving source ?
3966 17 : err += readVal( inrec, String("distance"), distance );
3967 17 : err += readVal( inrec, String("tracksource"), trackSource );
3968 17 : err += readVal( inrec, String("trackdir"), trackDir );
3969 :
3970 : // The extra params for WB-AWP
3971 17 : err += readVal( inrec, String("aterm"), aTermOn );
3972 17 : err += readVal( inrec, String("psterm"), psTermOn );
3973 17 : err += readVal( inrec, String("mterm"), mTermOn );
3974 17 : err += readVal( inrec, String("wbawp"), wbAWP );
3975 17 : err += readVal( inrec, String("cfcache"), cfCache );
3976 17 : err += readVal( inrec, String("usepointing"), usePointing );
3977 17 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3978 17 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3979 17 : err += readVal( inrec, String("conjbeams"), conjBeams );
3980 17 : err += readVal( inrec, String("computepastep"), computePAStep );
3981 17 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3982 :
3983 : // The extra params for single-dish
3984 17 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3985 17 : err += readVal( inrec, String("convertfirst"), convertFirst );
3986 17 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3987 17 : err += readVal( inrec, String("convsupport"), convSupport );
3988 17 : err += readVal( inrec, String("truncate"), truncateSize );
3989 17 : err += readVal( inrec, String("gwidth"), gwidth );
3990 17 : err += readVal( inrec, String("jwidth"), jwidth );
3991 17 : err += readVal( inrec, String("minweight"), minWeight );
3992 17 : 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 17 : if (ftmachine=="awprojectft" && cfCache=="") {
3998 0 : cfCache = imageName + ".cf";
3999 : }
4000 :
4001 17 : if ( ftmachine=="awprojectft" &&
4002 17 : 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 17 : err += verify();
4011 :
4012 17 : } catch(AipsError &x) {
4013 0 : err = err + x.getMesg() + "\n";
4014 0 : }
4015 :
4016 17 : if (err.length()>0) {
4017 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4018 : }
4019 :
4020 17 : }
4021 :
4022 17 : String SynthesisParamsGrid::verify() const
4023 : {
4024 17 : 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 17 : if ( imageName == "" ) {
4032 0 : err += "Please supply an image name\n";
4033 : }
4034 :
4035 0 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4036 17 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4037 17 : (ftmachine != "mawprojectft") &&
4038 0 : (ftmachine != "sd"))
4039 : {
4040 : err += "Invalid ftmachine name. Must be one of"
4041 : " 'gridft', 'wprojectft',"
4042 : " 'mosaicft', 'awprojectft',"
4043 : " 'mawpojectft', 'protoft',"
4044 0 : " 'sd'\n";
4045 : }
4046 :
4047 :
4048 34 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4049 17 : ( 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 17 : 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 17 : 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 17 : 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 17 : 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 17 : 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 17 : 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 17 : 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 17 : if ( ftmachine == "sd" ) {
4115 0 : if ( convertFirst != "always" and
4116 0 : convertFirst != "never" and
4117 0 : convertFirst != "auto" ) {
4118 0 : err += "convertfirst parameter: illegal value: '" + convertFirst + "'."
4119 0 : " Allowed values: 'always', 'never', 'auto'.\n";
4120 : }
4121 : }
4122 :
4123 17 : return err;
4124 0 : }
4125 :
4126 68 : void SynthesisParamsGrid::setDefaults()
4127 : {
4128 68 : imageName="";
4129 : // FTMachine parameters
4130 : //ftmachine="GridFT";
4131 68 : ftmachine="gridft";
4132 68 : gridder=ftmachine;
4133 68 : padding=1.2;
4134 68 : useAutoCorr=false;
4135 68 : useDoublePrec=true;
4136 68 : wprojplanes=1;
4137 68 : convFunc="SF";
4138 68 : vpTable="";
4139 :
4140 : // facets
4141 68 : facets=1;
4142 :
4143 : // chanchunks
4144 68 : chanchunks=1;
4145 :
4146 : // Spectral Axis interpolation
4147 68 : interpolation=String("nearest");
4148 :
4149 : //mosaic use pointing
4150 68 : usePointing=false;
4151 : // Moving phase center ?
4152 68 : distance=Quantity(0,"m");
4153 68 : trackSource=false;
4154 68 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4155 :
4156 : // The extra params for WB-AWP
4157 68 : aTermOn = true;
4158 68 : psTermOn = true;
4159 68 : mTermOn = false;
4160 68 : wbAWP = true;
4161 68 : cfCache = "";
4162 68 : usePointing = false;
4163 68 : pointingOffsetSigDev.resize(0);
4164 : // pointingOffsetSigDev.set(30.0);
4165 68 : doPBCorr = true;
4166 68 : conjBeams = true;
4167 68 : computePAStep=360.0;
4168 68 : rotatePAStep=5.0;
4169 :
4170 : // extra params for single-dish
4171 68 : pointingDirCol = "";
4172 68 : convertFirst = "never";
4173 68 : skyPosThreshold = 0.0;
4174 68 : convSupport = -1;
4175 68 : truncateSize = Quantity(-1.0);
4176 68 : gwidth = Quantity(-1.0);
4177 68 : jwidth = Quantity(-1.0);
4178 68 : minWeight = 0.0;
4179 68 : clipMinMax = False;
4180 :
4181 : // Mapper type
4182 68 : mType = String("default");
4183 :
4184 68 : }
4185 :
4186 0 : Record SynthesisParamsGrid::toRecord() const
4187 : {
4188 0 : Record gridpar;
4189 :
4190 0 : gridpar.define("imagename", imageName);
4191 : // FTMachine params
4192 0 : gridpar.define("padding", padding);
4193 0 : gridpar.define("useautocorr",useAutoCorr );
4194 0 : gridpar.define("usedoubleprec", useDoublePrec);
4195 0 : gridpar.define("wprojplanes", wprojplanes);
4196 0 : gridpar.define("convfunc", convFunc);
4197 0 : gridpar.define("vptable", vpTable);
4198 :
4199 0 : gridpar.define("facets", facets);
4200 0 : gridpar.define("chanchunks", chanchunks);
4201 :
4202 0 : gridpar.define("interpolation",interpolation);
4203 :
4204 0 : gridpar.define("distance", QuantityToString(distance));
4205 0 : gridpar.define("tracksource", trackSource);
4206 0 : gridpar.define("trackdir", MDirectionToString( trackDir ));
4207 :
4208 0 : gridpar.define("aterm",aTermOn );
4209 0 : gridpar.define("psterm",psTermOn );
4210 0 : gridpar.define("mterm",mTermOn );
4211 0 : gridpar.define("wbawp", wbAWP);
4212 0 : gridpar.define("cfcache", cfCache);
4213 0 : gridpar.define("usepointing",usePointing );
4214 0 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
4215 0 : gridpar.define("dopbcorr", doPBCorr);
4216 0 : gridpar.define("conjbeams",conjBeams );
4217 0 : gridpar.define("computepastep", computePAStep);
4218 0 : gridpar.define("rotatepastep", rotatePAStep);
4219 :
4220 0 : gridpar.define("pointingcolumntouse", pointingDirCol );
4221 0 : gridpar.define("convertfirst", convertFirst );
4222 0 : gridpar.define("skyposthreshold", skyPosThreshold );
4223 0 : gridpar.define("convsupport", convSupport );
4224 0 : gridpar.define("truncate", QuantityToString(truncateSize) );
4225 0 : gridpar.define("gwidth", QuantityToString(gwidth) );
4226 0 : gridpar.define("jwidth", QuantityToString(jwidth) );
4227 0 : gridpar.define("minweight", minWeight );
4228 0 : gridpar.define("clipminmax", clipMinMax );
4229 :
4230 0 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
4231 : /// else gridpar.define("deconvolver","singleterm");
4232 :
4233 0 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
4234 0 : else gridpar.define("gridder", gridder);
4235 :
4236 : // gridpar.define("mtype", mType);
4237 :
4238 0 : return gridpar;
4239 0 : }
4240 :
4241 :
4242 :
4243 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
4244 :
4245 : /////////////////////// Deconvolver Parameters
4246 :
4247 34 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4248 : {
4249 34 : setDefaults();
4250 34 : }
4251 :
4252 34 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4253 : {
4254 34 : }
4255 :
4256 :
4257 34 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4258 : {
4259 34 : setDefaults();
4260 :
4261 34 : String err("");
4262 :
4263 : try
4264 : {
4265 :
4266 34 : err += readVal( inrec, String("imagename"), imageName );
4267 34 : 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 34 : if( inrec.isDefined("startmodel") )
4273 : {
4274 34 : if( inrec.dataType("startmodel")==TpString )
4275 : {
4276 17 : String onemodel;
4277 17 : err += readVal( inrec, String("startmodel"), onemodel );
4278 17 : if( onemodel.length()>0 )
4279 : {
4280 0 : startModel.resize(1);
4281 0 : startModel[0] = onemodel;
4282 : }
4283 17 : else {startModel.resize();}
4284 17 : }
4285 34 : else if( inrec.dataType("startmodel")==TpArrayString ||
4286 17 : inrec.dataType("startmodel")==TpArrayBool)
4287 : {
4288 17 : 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 34 : err += readVal( inrec, String("id"), deconvolverId );
4297 34 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4298 :
4299 34 : err += readVal( inrec, String("scales"), scales );
4300 34 : err += readVal( inrec, String("scalebias"), scalebias );
4301 :
4302 34 : err += readVal( inrec, String("usemask"), maskType );
4303 34 : if( maskType=="auto-thresh" )
4304 : {
4305 0 : autoMaskAlgorithm = "thresh";
4306 : }
4307 34 : else if( maskType=="auto-thesh2" )
4308 : {
4309 0 : autoMaskAlgorithm = "thresh2";
4310 : }
4311 34 : else if( maskType=="auto-onebox" )
4312 : {
4313 0 : autoMaskAlgorithm = "onebox";
4314 : }
4315 34 : else if( maskType=="user" || maskType=="pb" )
4316 : {
4317 34 : autoMaskAlgorithm = "";
4318 : }
4319 :
4320 :
4321 34 : if( inrec.isDefined("mask") )
4322 : {
4323 34 : if( inrec.dataType("mask")==TpString )
4324 : {
4325 34 : err+= readVal( inrec, String("mask"), maskString );
4326 : }
4327 0 : else if( inrec.dataType("mask")==TpArrayString )
4328 : {
4329 0 : err+= readVal( inrec, String("mask"), maskList );
4330 : }
4331 : }
4332 :
4333 34 : if( inrec.isDefined("pbmask") )
4334 : {
4335 34 : err += readVal( inrec, String("pbmask"), pbMask );
4336 : }
4337 34 : if( inrec.isDefined("maskthreshold") )
4338 : {
4339 34 : if( inrec.dataType("maskthreshold")==TpString )
4340 : {
4341 34 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4342 : //deal with the case a string is a float value without unit
4343 34 : Quantity testThresholdString;
4344 34 : Quantity::read(testThresholdString,maskThreshold);
4345 34 : if( testThresholdString.getUnit()=="" )
4346 : {
4347 34 : if(testThresholdString.getValue()<1.0)
4348 : {
4349 34 : fracOfPeak = testThresholdString.getValue();
4350 34 : maskThreshold=String("");
4351 : }
4352 : }
4353 34 : }
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 34 : if( inrec.isDefined("maskresolution") )
4372 : {
4373 34 : if( inrec.dataType("maskresolution")==TpString )
4374 : {
4375 34 : err += readVal(inrec, String("maskresolution"), maskResolution );
4376 : //deal with the case a string is a float value without unit
4377 34 : Quantity testResolutionString;
4378 34 : Quantity::read(testResolutionString,maskResolution);
4379 34 : if( testResolutionString.getUnit()=="" )
4380 : {
4381 34 : maskResByBeam = testResolutionString.getValue();
4382 34 : maskResolution=String("");
4383 : }
4384 34 : }
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 34 : if( inrec.isDefined("nmask") )
4398 : {
4399 34 : if( inrec.dataType("nmask")==TpInt )
4400 : {
4401 34 : err+= readVal(inrec, String("nmask"), nMask );
4402 : }
4403 : else
4404 : {
4405 0 : err+= "nmask must be an integer\n";
4406 : }
4407 : }
4408 34 : if( inrec.isDefined("autoadjust") )
4409 : {
4410 17 : if( inrec.dataType("autoadjust")==TpBool )
4411 : {
4412 17 : 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 34 : if (inrec.isDefined("fusedthreshold"))
4421 : {
4422 34 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4423 34 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4424 : else
4425 0 : err += "fusedthreshold must be a float or double";
4426 : }
4427 34 : if (inrec.isDefined("specmode"))
4428 : {
4429 34 : if(inrec.dataType("specmode") == TpString)
4430 34 : 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 34 : if (inrec.isDefined("largestscale"))
4436 : {
4437 34 : if (inrec.dataType("largestscale") == TpInt)
4438 34 : 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 34 : if( inrec.isDefined("sidelobethreshold"))
4444 : {
4445 34 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4446 : {
4447 34 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4448 : }
4449 : else
4450 : {
4451 0 : err+= "sidelobethreshold must be a float or double";
4452 : }
4453 : }
4454 :
4455 34 : if( inrec.isDefined("noisethreshold"))
4456 : {
4457 34 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4458 : {
4459 34 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4460 : }
4461 : else
4462 : {
4463 0 : err+= "noisethreshold must be a float or double";
4464 : }
4465 : }
4466 34 : if( inrec.isDefined("lownoisethreshold"))
4467 : {
4468 34 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4469 : {
4470 34 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4471 : }
4472 : else
4473 : {
4474 0 : err+= "lownoisethreshold must be a float or double";
4475 : }
4476 : }
4477 34 : if( inrec.isDefined("negativethreshold"))
4478 : {
4479 34 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4480 : {
4481 34 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4482 : }
4483 : else
4484 : {
4485 0 : err+= "negativethreshold must be a float or double";
4486 : }
4487 : }
4488 34 : if( inrec.isDefined("smoothfactor"))
4489 : {
4490 34 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4491 : {
4492 34 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4493 : }
4494 : else
4495 : {
4496 0 : err+= "smoothfactor must be a float or double";
4497 : }
4498 : }
4499 34 : if( inrec.isDefined("minbeamfrac"))
4500 : {
4501 34 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4502 : {
4503 34 : 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 34 : if( inrec.isDefined("cutthreshold"))
4517 : {
4518 34 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4519 : {
4520 34 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4521 : }
4522 : else {
4523 0 : err+= "cutthreshold must be a float or double";
4524 : }
4525 : }
4526 34 : if( inrec.isDefined("growiterations"))
4527 : {
4528 34 : if (inrec.dataType("growiterations")==TpInt) {
4529 34 : err+= readVal(inrec, String("growiterations"), growIterations );
4530 : }
4531 : else {
4532 0 : err+= "growiterations must be an integer\n";
4533 : }
4534 : }
4535 34 : if( inrec.isDefined("dogrowprune"))
4536 : {
4537 34 : if (inrec.dataType("dogrowprune")==TpBool) {
4538 34 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4539 : }
4540 : else {
4541 0 : err+= "dogrowprune must be a bool\n";
4542 : }
4543 : }
4544 34 : if( inrec.isDefined("minpercentchange"))
4545 : {
4546 34 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4547 34 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4548 : }
4549 : else {
4550 0 : err+= "minpercentchange must be a float or double";
4551 : }
4552 : }
4553 34 : if( inrec.isDefined("verbose"))
4554 : {
4555 34 : if (inrec.dataType("verbose")==TpBool ) {
4556 34 : err+= readVal(inrec, String("verbose"), verbose);
4557 : }
4558 : else {
4559 0 : err+= "verbose must be a bool";
4560 : }
4561 : }
4562 34 : if( inrec.isDefined("fastnoise"))
4563 : {
4564 34 : if (inrec.dataType("fastnoise")==TpBool ) {
4565 34 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4566 : }
4567 : else {
4568 0 : err+= "fastnoise must be a bool";
4569 : }
4570 : }
4571 34 : if( inrec.isDefined("nsigma") )
4572 : {
4573 34 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4574 34 : 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 34 : if( inrec.isDefined("noRequireSumwt") )
4587 : {
4588 17 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4589 17 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4590 : }
4591 : else {
4592 0 : err+= "noRequireSumwt must be a bool";
4593 : }
4594 : }
4595 34 : if( inrec.isDefined("fullsummary") )
4596 : {
4597 34 : if (inrec.dataType("fullsummary")==TpBool) {
4598 34 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4599 : }
4600 : else {
4601 0 : err+= "fullsummary must be a bool";
4602 : }
4603 : }
4604 34 : if( inrec.isDefined("restoringbeam") )
4605 : {
4606 17 : String errinfo("");
4607 : try {
4608 :
4609 17 : 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 17 : 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 17 : }// if isdefined(restoringbeam)
4658 :
4659 34 : if( inrec.isDefined("interactive") )
4660 : {
4661 34 : if( inrec.dataType("interactive")==TpBool )
4662 34 : {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 34 : err += verify();
4670 :
4671 : }
4672 0 : catch(AipsError &x)
4673 : {
4674 0 : err = err + x.getMesg() + "\n";
4675 0 : }
4676 :
4677 34 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4678 :
4679 34 : }
4680 :
4681 34 : String SynthesisParamsDeconv::verify() const
4682 : {
4683 34 : String err;
4684 :
4685 34 : 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 34 : if( maskString.length()>0 )
4689 : {
4690 34 : File fp( imageName+".mask" );
4691 34 : 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 34 : }
4693 :
4694 34 : if( pbMask >= 1.0)
4695 0 : {err += "pbmask must be < 1.0 \n"; }
4696 34 : else if( pbMask < 0.0)
4697 0 : {err += "pbmask must be a positive value \n"; }
4698 :
4699 34 : 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 34 : if ( fracOfPeak >= 1.0)
4708 0 : {err += "fracofpeak must be < 1.0 \n"; }
4709 34 : else if ( fracOfPeak < 0.0)
4710 0 : {err += "fracofpeak must be a positive value \n"; }
4711 :
4712 34 : return err;
4713 0 : }
4714 :
4715 68 : void SynthesisParamsDeconv::setDefaults()
4716 : {
4717 68 : imageName="";
4718 68 : algorithm="hogbom";
4719 68 : startModel=Vector<String>(0);
4720 68 : deconvolverId=0;
4721 68 : nTaylorTerms=1;
4722 68 : scales.resize(1); scales[0]=0.0;
4723 68 : scalebias=0.6;
4724 68 : maskType="none";
4725 68 : maskString="";
4726 68 : maskList.resize(1); maskList[0]="";
4727 68 : pbMask=0.0;
4728 68 : autoMaskAlgorithm="thresh";
4729 68 : maskThreshold="";
4730 68 : maskResolution="";
4731 68 : fracOfPeak=0.0;
4732 68 : nMask=0;
4733 68 : interactive=false;
4734 68 : autoAdjust=False;
4735 68 : fusedThreshold = 0.0;
4736 68 : specmode="mfs";
4737 68 : largestscale = -1;
4738 68 : }
4739 :
4740 17 : Record SynthesisParamsDeconv::toRecord() const
4741 : {
4742 17 : Record decpar;
4743 :
4744 17 : decpar.define("imagename", imageName);
4745 17 : decpar.define("deconvolver", algorithm);
4746 17 : decpar.define("startmodel",startModel);
4747 17 : decpar.define("id",deconvolverId);
4748 17 : decpar.define("nterms",nTaylorTerms);
4749 17 : decpar.define("scales",scales);
4750 17 : decpar.define("scalebias",scalebias);
4751 17 : decpar.define("usemask",maskType);
4752 17 : decpar.define("fusedthreshold", fusedThreshold);
4753 17 : decpar.define("specmode", specmode);
4754 17 : decpar.define("largestscale", largestscale);
4755 17 : if( maskList.nelements()==1 && maskList[0]=="")
4756 : {
4757 17 : decpar.define("mask",maskString);
4758 : }
4759 : else {
4760 0 : decpar.define("mask",maskList);
4761 : }
4762 17 : decpar.define("pbmask",pbMask);
4763 17 : if (fracOfPeak > 0.0)
4764 : {
4765 0 : decpar.define("maskthreshold",fracOfPeak);
4766 : }
4767 : else
4768 : {
4769 17 : decpar.define("maskthreshold",maskThreshold);
4770 : }
4771 17 : decpar.define("maskresolution",maskResolution);
4772 17 : decpar.define("nmask",nMask);
4773 17 : decpar.define("autoadjust",autoAdjust);
4774 17 : decpar.define("sidelobethreshold",sidelobeThreshold);
4775 17 : decpar.define("noisethreshold",noiseThreshold);
4776 17 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4777 17 : decpar.define("negativethreshold",negativeThreshold);
4778 17 : decpar.define("smoothfactor",smoothFactor);
4779 17 : decpar.define("minbeamfrac",minBeamFrac);
4780 17 : decpar.define("cutthreshold",cutThreshold);
4781 17 : decpar.define("growiterations",growIterations);
4782 17 : decpar.define("dogrowprune",doGrowPrune);
4783 17 : decpar.define("minpercentchange",minPercentChange);
4784 17 : decpar.define("verbose", verbose);
4785 17 : decpar.define("fastnoise", fastnoise);
4786 17 : decpar.define("interactive",interactive);
4787 17 : decpar.define("nsigma",nsigma);
4788 17 : decpar.define("noRequireSumwt",noRequireSumwt);
4789 17 : decpar.define("fullsummary",fullsummary);
4790 :
4791 17 : return decpar;
4792 0 : }
4793 :
4794 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4795 :
4796 :
4797 : } //# NAMESPACE CASA - END
4798 :
|
