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 1446 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 1446 : }
89 :
90 1446 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 1446 : }
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 514533 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 514533 : Int N=vb.nRows(),M=-1;
108 2269669 : for(Int i=0;i<N;i++)
109 : {
110 2269669 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 514533 : {M++;break;}
112 : }
113 514533 : 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 3698 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 3698 : Int n=npix;
121 :
122 3698 : if( n%2 !=0 ){ n+= 1; }
123 :
124 3698 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 24001 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 20303 : if( fac[k]>5 )
129 : {
130 147 : val = fac[k];
131 301 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 147 : fac[k] = val;
133 : }
134 : }
135 3698 : newlarge=product(fac);
136 3955 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 272 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 3683 : return newlarge;
141 3698 : }
142 :
143 : // Return the prime factors of the given number
144 4271 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 4271 : Vector<uInt> factors;
147 :
148 4271 : Int lastresult = n;
149 4271 : Int sqlast = int(sqrt(n))+1;
150 :
151 4271 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 4271 : Int c=2;
153 : while(1)
154 : {
155 25980 : if( lastresult == 1 || c > sqlast ) { break; }
156 21709 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 31351 : if( c>sqlast){ c=lastresult; break; }
160 28082 : if( lastresult % c == 0 ) { break; }
161 9642 : c += 1;
162 : }
163 21709 : factors.resize( factors.nelements()+1, true );
164 21709 : factors[ factors.nelements()-1 ] = c;
165 21709 : lastresult /= c;
166 : }
167 4271 : 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 4271 : return factors;
185 0 : }
186 :
187 :
188 35 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 70 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 35 : 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 35 : std::shared_ptr<SIImageStore> imstore;
199 35 : if( nterms>1 )
200 10 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 25 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 35 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 35 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 35 : imstore->printBeamSet();
210 :
211 35 : imstore->releaseLocks();
212 :
213 35 : return true;
214 35 : }
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 18 : Bool SynthesisUtilMethods::taylorCoeffsToCube(const String& cubename,const String& mtname, const Int nterms, const String& reffreq)
222 : {
223 36 : LogIO os(LogOrigin("SynthesisUtilMethods", "taylorCoeffsToCube"));
224 :
225 : // Set up imstores
226 18 : CountedPtr<SIImageStore> cube_imstore;
227 18 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
228 :
229 18 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
230 18 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
231 :
232 : // Check that .model exists.
233 : try{
234 18 : cube_imstore->model();
235 54 : for(Int i=0;i<nterms;i++)
236 36 : 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 18 : IPosition cube_shp( cube_imstore->model()->shape() );
245 18 : IPosition mt_shp( mt_imstore->model(0)->shape() );
246 18 : 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 18 : Quantity reffreq_qa;
252 18 : Quantity::read( reffreq_qa, reffreq );
253 18 : Double refval = reffreq_qa.getValue("Hz");
254 :
255 : //cout << "ref freq : " << refval << endl;
256 :
257 : //Get the frequency list for the cube
258 18 : CoordinateSystem csys ( cube_imstore->getCSys() );
259 18 : Vector<Double> freqlist( cube_shp[3] );
260 :
261 72 : for(uInt i=0; i<csys.nCoordinates(); i++)
262 : {
263 54 : if( csys.type(i) == Coordinate::SPECTRAL )
264 : {
265 18 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
266 :
267 108 : for(Int ch=0;ch<cube_shp[3];ch++)
268 : {
269 : Double freq;
270 90 : Bool ret = speccoord.toWorld( freq, ch );
271 90 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
272 90 : freqlist[ch] = freq;
273 :
274 : //cout << "freq " << ch << " is " << freq << endl;
275 : }
276 :
277 18 : }
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 36 : for(Int pol=0; pol<cube_shp[2]; pol++)
286 : {
287 18 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(nterms);
288 54 : for(Int i=0;i<nterms;i++)
289 : {
290 72 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
291 36 : 0, cube_shp[3],
292 36 : pol, cube_shp[2],
293 108 : *mt_imstore->model(i) );
294 : }
295 :
296 108 : for(Int chan=0; chan<cube_shp[3]; chan++)
297 : {
298 180 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
299 90 : chan, cube_shp[3],
300 90 : pol, cube_shp[2],
301 180 : *cube_imstore->model() );
302 :
303 90 : Double wt = (freqlist[chan] - refval) / refval;
304 :
305 90 : cube_subim->set(0.0);
306 270 : for(Int tt=0;tt<nterms;tt++)
307 : {
308 180 : Double fac = pow(wt,tt);
309 360 : LatticeExpr<Float> oneterm = LatticeExpr<Float>( *cube_subim + (fac) * (*mt_subims[tt])) ;
310 180 : cube_subim->copyData(oneterm);
311 180 : }
312 :
313 :
314 90 : }//for chan
315 18 : }// for pol
316 :
317 18 : return True;
318 :
319 18 : }//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 53 : Bool SynthesisUtilMethods::cubeToTaylorSum(const String& cubename,const String& mtname, const Int nterms, const String& reffreq, const Int imtype, const Float pblimit)
330 : {
331 106 : LogIO os(LogOrigin("SynthesisUtilMethods", "cubeToTaylorSum"));
332 :
333 : //cout << "imtype : " << imtype << endl;
334 53 : 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 53 : CountedPtr<SIImageStore> cube_imstore;
341 53 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
342 :
343 53 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
344 53 : 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 53 : Float maxPB=1.0;
348 53 : if(imtype < 2){
349 70 : LatticeExprNode elnod( max( *(mt_imstore->pb(0)) ) );
350 35 : maxPB=elnod.getFloat();
351 35 : if(maxPB == 0.0){
352 0 : throw(AipsError("Programmers error: should do tt psf images after making average PB"));
353 :
354 : }
355 :
356 35 : }
357 : // cerr << "imtype " << imtype << " MAX PB " << maxPB << endl;
358 : // If dopsf=True, calculate 2n-1 terms.
359 53 : Int out_nterms=nterms; // for residual
360 53 : 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 53 : 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 53 : 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 53 : switch(imtype)
367 : {
368 9 : case 0: use_cube=cube_imstore->psf();break;
369 26 : case 1: use_cube=cube_imstore->residual();break;
370 9 : case 2: use_cube=cube_imstore->pb();break;
371 9 : case 3: use_cube=cube_imstore->sumwt();break;
372 : }
373 53 : cube_imstore->sumwt();
374 168 : for(Int i=0;i<out_nterms;i++)
375 : {
376 115 : switch(imtype)
377 : {
378 27 : case 0:mt_imstore->psf(i);break;
379 52 : case 1:mt_imstore->residual(i);break;
380 9 : case 2:mt_imstore->pb(i);break;
381 27 : 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 106 : IPosition cube_shp( cube_imstore->psf()->shape() );
392 106 : IPosition mt_shp( mt_imstore->psf(0)->shape() );
393 53 : 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 106 : Quantity reffreq_qa;
399 53 : Quantity::read( reffreq_qa, reffreq );
400 53 : Double refval = reffreq_qa.getValue("Hz");
401 :
402 : //cout << "ref freq : " << refval << endl;
403 :
404 : //Get the frequency list for the cube
405 106 : CoordinateSystem csys ( cube_imstore->getCSys() );
406 106 : Vector<Double> freqlist( cube_shp[3] );
407 :
408 212 : for(uInt i=0; i<csys.nCoordinates(); i++)
409 : {
410 159 : if( csys.type(i) == Coordinate::SPECTRAL )
411 : {
412 53 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
413 :
414 324 : for(Int ch=0;ch<cube_shp[3];ch++)
415 : {
416 : Double freq;
417 271 : Bool ret = speccoord.toWorld( freq, ch );
418 271 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
419 271 : freqlist[ch] = freq;
420 : // cout << "freq " << ch << " is " << freq << endl;
421 : }
422 :
423 53 : }
424 : }
425 :
426 :
427 : // Reset the Taylor Sum values to zero.
428 168 : for(Int i=0;i<out_nterms;i++)
429 : {
430 115 : switch(imtype)
431 : {
432 27 : case 0:mt_imstore->psf(i)->set(0.0);break;
433 52 : case 1:mt_imstore->residual(i)->set(0.0);break;
434 9 : case 2:mt_imstore->pb(i)->set(0.0);break;
435 27 : case 3:mt_imstore->sumwt(i)->set(0.0);break;
436 : }
437 : }
438 :
439 : // Get the sumwt spectrum.
440 106 : Array<Float> lsumwt;
441 53 : cube_imstore->sumwt()->get(lsumwt, False);
442 :
443 : // Sum the weights ( or just use accumulate...)
444 106 : LatticeExprNode msum( sum( *cube_imstore->sumwt() ) );
445 53 : 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 106 : for(Int pol=0; pol<cube_shp[2]; pol++)
453 : {
454 53 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(out_nterms);
455 168 : for(Int i=0;i<out_nterms;i++)
456 : {
457 115 : switch(imtype)
458 : {
459 27 : case 0:use_mt=mt_imstore->psf(i);break;
460 52 : case 1:use_mt=mt_imstore->residual(i);break;
461 9 : case 2:use_mt=mt_imstore->pb(i);break;
462 27 : case 3:use_mt=mt_imstore->sumwt(i);break;
463 : }
464 345 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
465 115 : 0, cube_shp[3],
466 115 : pol, cube_shp[2],
467 230 : *use_mt );
468 : }
469 :
470 324 : for(Int chan=0; chan<cube_shp[3]; chan++)
471 : {
472 271 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
473 271 : chan, cube_shp[3],
474 271 : pol, cube_shp[2],
475 271 : *use_cube);
476 271 : if(imtype < 2){
477 358 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
478 179 : chan, cube_shp[3],
479 179 : pol, cube_shp[2],
480 358 : *(cube_imstore->pb()));
481 358 : CountedPtr<ImageInterface<Float> > tmplat = new TempImage<Float>(cube_subim->shape(), cube_subim->coordinates());
482 179 : tmplat->copyData(LatticeExpr<Float>((*pb_subim) *(*cube_subim)));
483 179 : cube_subim = tmplat;
484 :
485 179 : }
486 :
487 271 : IPosition pos(4,0,0,pol,chan);
488 :
489 271 : Double wt = (freqlist[chan] - refval) / refval;
490 :
491 859 : for(Int tt=0;tt<out_nterms;tt++)
492 : {
493 588 : Double fac = pow(wt,tt);
494 : //cerr << "BEF accum " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
495 1176 : LatticeExpr<Float> eachterm = LatticeExpr<Float>( (*mt_subims[tt]) + ((fac) * (*cube_subim) * lsumwt(pos))) ;
496 588 : 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 588 : }
499 :
500 :
501 271 : }//for chan
502 :
503 :
504 : // Divide by sum of weights.
505 168 : 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 115 : LatticeExpr<Float> eachterm;
510 115 : if (imtype < 2) {
511 79 : eachterm = LatticeExpr<Float>( iif( (*(mt_imstore->pb(0))) > pblimit , (*mt_subims[tt]) / wtsum/(*(mt_imstore->pb(0))), 0.0));
512 : }
513 : else{
514 36 : eachterm = LatticeExpr<Float>( (*mt_subims[tt]) / wtsum ) ;
515 : }
516 115 : mt_subims[tt]->copyData(eachterm);
517 :
518 115 : mt_subims[tt]->flush();
519 : // cerr << "aft div : " << max(mt_subims[tt]->get()) << endl;
520 115 : }
521 :
522 53 : }// 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 53 : if( imtype==2 )
527 : {
528 9 : mt_imstore->removeMask( mt_imstore->pb(0) );
529 : {
530 : //MSK//
531 18 : LatticeExpr<Bool> pbmask( iif( *mt_imstore->pb(0) > fabs(pblimit) , True , False ) );
532 : //MSK//
533 9 : mt_imstore->createMask( pbmask, mt_imstore->pb(0) );
534 9 : mt_imstore->pb(0)->pixelMask().unlock();
535 9 : }
536 :
537 : }
538 :
539 106 : return True;
540 :
541 53 : }//end of func
542 :
543 :
544 26 : Bool SynthesisUtilMethods::removeFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
545 : {
546 52 : LogIO os(LogOrigin("SynthesisUtilMethods", "removeFreqDepPB"));
547 :
548 : // Set up imstores
549 26 : CountedPtr<SIImageStore> cube_imstore;
550 26 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
551 :
552 26 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
553 26 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
554 :
555 : try{
556 26 : cube_imstore->residual(); // Residual Cube
557 26 : cube_imstore->pb(); // PB cube
558 26 : 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 26 : IPosition cube_shp( cube_imstore->residual()->shape() );
567 26 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
568 26 : 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 52 : for(Int pol=0; pol<cube_shp[2]; pol++)
574 : {
575 :
576 52 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
577 26 : 0, cube_shp[3],
578 26 : pol, cube_shp[2],
579 52 : (*mt_imstore->pb(0)) );
580 :
581 26 : LatticeExprNode mtpbmax( max( *mt_subim ) );
582 26 : Float mtpbmaxval = mtpbmax.getFloat();
583 26 : if(mtpbmaxval <=0.0){os << LogIO::WARN << "pb.tt0 max is < or = zero" << LogIO::POST;}
584 :
585 :
586 159 : for(Int chan=0; chan<cube_shp[3]; chan++)
587 : {
588 266 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
589 133 : chan, cube_shp[3],
590 133 : pol, cube_shp[2],
591 266 : *cube_imstore->residual() );
592 266 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
593 133 : chan, cube_shp[3],
594 133 : pol, cube_shp[2],
595 266 : *cube_imstore->pb() );
596 :
597 133 : LatticeExprNode pbmax( max( *pb_subim ) );
598 133 : Float pbmaxval = pbmax.getFloat();
599 133 : if( pbmaxval<=0.0 )
600 : {
601 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
602 : }
603 : else
604 : {
605 266 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*mt_subim)*(*(cube_subim))/(*(pb_subim)) , 0.0 ) );
606 133 : cube_subim->copyData( thepbcor );
607 133 : }// if not zero
608 :
609 133 : }//for chan
610 26 : }// for pol
611 :
612 26 : return True;
613 :
614 26 : }//end of func
615 :
616 :
617 :
618 18 : Bool SynthesisUtilMethods::applyFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
619 : {
620 36 : LogIO os(LogOrigin("SynthesisUtilMethods", "applyFreqDepPB"));
621 :
622 : // Set up imstores
623 18 : CountedPtr<SIImageStore> cube_imstore;
624 18 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
625 :
626 18 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
627 18 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
628 :
629 : try{
630 18 : cube_imstore->model(); // Model Cube
631 18 : cube_imstore->pb(); // PB cube
632 18 : 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 18 : IPosition cube_shp( cube_imstore->model()->shape() );
641 18 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
642 18 : 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 36 : for(Int pol=0; pol<cube_shp[2]; pol++)
648 : {
649 36 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
650 18 : 0, cube_shp[3],
651 18 : pol, cube_shp[2],
652 36 : (*mt_imstore->pb(0)) );
653 :
654 18 : LatticeExprNode mtpbmax( max( *mt_subim ) );
655 18 : Float mtpbmaxval = mtpbmax.getFloat();
656 18 : 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 108 : for(Int chan=0; chan<cube_shp[3]; chan++)
662 : {
663 180 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
664 90 : chan, cube_shp[3],
665 90 : pol, cube_shp[2],
666 180 : *cube_imstore->model() );
667 180 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
668 90 : chan, cube_shp[3],
669 90 : pol, cube_shp[2],
670 180 : *cube_imstore->pb() );
671 :
672 :
673 90 : LatticeExprNode pbmax( max( *pb_subim ) );
674 90 : Float pbmaxval = pbmax.getFloat();
675 90 : if( pbmaxval<=0.0 )
676 : {
677 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
678 : }
679 : else
680 : {
681 180 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*(cube_subim)) *(*(pb_subim)) / (*mt_subim) , 0.0 ) );
682 90 : cube_subim->copyData( thepbcor );
683 90 : }// if not zero
684 :
685 90 : }//for chan
686 : }// if mtpb >0
687 18 : }// for pol
688 :
689 18 : return True;
690 :
691 18 : }//end of func
692 :
693 :
694 :
695 :
696 : /***make a record of synthesisimager::weight parameters***/
697 1630 : 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 1630 : Record outRec;
705 1630 : outRec.define("type", type);
706 1630 : outRec.define("rmode", rmode);
707 1630 : Record quantRec;
708 1630 : QuantumHolder(noise).toRecord(quantRec);
709 1630 : outRec.defineRecord("noise", quantRec);
710 1630 : outRec.define("robust", robust);
711 1630 : QuantumHolder(fieldofview).toRecord(quantRec);
712 1630 : outRec.defineRecord("fieldofview", quantRec);
713 1630 : outRec.define("npixels", npixels);
714 1630 : outRec.define("multifield", multiField);
715 1630 : outRec.define("usecubebriggs", useCubeBriggs);
716 1630 : outRec.define("filtertype", filtertype);
717 1630 : QuantumHolder(filterbmaj).toRecord(quantRec);
718 1630 : outRec.defineRecord("filterbmaj", quantRec);
719 1630 : QuantumHolder(filterbmin).toRecord(quantRec);
720 1630 : outRec.defineRecord("filterbmin", quantRec);
721 1630 : QuantumHolder(filterbpa).toRecord(quantRec);
722 1630 : outRec.defineRecord("filterbpa", quantRec);
723 1630 : outRec.define("fracBW", fracBW);
724 :
725 3260 : return outRec;
726 1630 : }
727 869 : 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 869 : QuantumHolder qh;
734 869 : String err;
735 869 : if(!inRec.isDefined("type"))
736 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
737 869 : inRec.get("type", type);
738 869 : inRec.get("rmode", rmode);
739 869 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
740 0 : throw(AipsError("Error in reading noise param"));
741 869 : noise=qh.asQuantity();
742 869 : inRec.get("robust", robust);
743 869 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
744 0 : throw(AipsError("Error in reading fieldofview param"));
745 869 : fieldofview=qh.asQuantity();
746 869 : inRec.get("npixels", npixels);
747 869 : inRec.get("multifield", multiField);
748 869 : inRec.get("usecubebriggs", useCubeBriggs);
749 869 : inRec.get("filtertype", filtertype);
750 869 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
751 0 : throw(AipsError("Error in reading filterbmaj param"));
752 869 : filterbmaj=qh.asQuantity();
753 869 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
754 0 : throw(AipsError("Error in reading filterbmin param"));
755 869 : filterbmin=qh.asQuantity();
756 869 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
757 0 : throw(AipsError("Error in reading filterbpa param"));
758 869 : filterbpa=qh.asQuantity();
759 869 : inRec.get("fracBW", fracBW);
760 :
761 :
762 :
763 869 : }
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 13 : 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 26 : LogIO os(LogOrigin("SynthesisUtilMethods", "adviseChanSel"));
818 13 : if(ms==String("")){
819 0 : throw(AipsError("Need a valid MS"));
820 : }
821 13 : spw.resize();
822 13 : start.resize();
823 13 : nchan.resize();
824 : try {
825 13 : 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 13 : MeasurementSet elms(ms, TableLock(TableLock::AutoNoReadLocking), Table::Old);
857 13 : MSSelection thisSelection;
858 13 : String spsel=spwselection;
859 13 : if(spsel=="")spsel="*";
860 13 : thisSelection.setSpwExpr(spsel);
861 13 : TableExprNode exprNode=thisSelection.toTableExprNode(&elms);
862 13 : Matrix<Int> chanlist=thisSelection.getChanList();
863 13 : if(chanlist.ncolumn() <3){
864 0 : freqStart=-1.0;
865 0 : freqEnd=-1.0;
866 0 : return false;
867 : }
868 13 : Vector<Int> elspw=chanlist.column(0);
869 13 : Vector<Int> elstart=chanlist.column(1);
870 26 : Vector<Int> elnchan=Vector<Int> (chanlist.column(2)-elstart)+1;
871 13 : 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 13 : MSUtil::getFreqRangeInSpw(freqStart, freqEnd, elspw, elstart, elnchan, elms, freqframe, field_id);
893 13 : }
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 13 : return true;
913 :
914 13 : }
915 :
916 26901 : 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 26901 : Bool isOn = false;
924 26901 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
925 26901 : if (!isOn)
926 26901 : 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 155 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
1398 : {
1399 155 : MVAngle mvRa=direction.getAngle().getValue()(0);
1400 155 : MVAngle mvDec=direction.getAngle().getValue()(1);
1401 155 : ostringstream oos;
1402 155 : oos << " ";
1403 155 : Int widthRA=20;
1404 155 : Int widthDec=20;
1405 155 : oos.setf(ios::left, ios::adjustfield);
1406 155 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
1407 155 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
1408 155 : oos << " "
1409 155 : << MDirection::showType(direction.getRefPtr()->getType());
1410 310 : return String(oos);
1411 155 : }
1412 :
1413 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1414 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1415 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
1416 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1417 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1418 :
1419 : // Read String from Record
1420 99791 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1421 : {
1422 99791 : if( rec.isDefined( id ) )
1423 : {
1424 92290 : String inval("");
1425 92290 : if( rec.dataType( id )==TpString )
1426 92290 : { 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 92290 : if(inval.length()>0){val=inval;}
1432 92290 : return String("");
1433 : }
1434 0 : else { return String(id + " must be a string\n"); }
1435 92290 : }
1436 7501 : else { return String("");}
1437 : }
1438 :
1439 : // Read Integer from Record
1440 25412 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1441 : {
1442 25412 : if( rec.isDefined( id ) )
1443 : {
1444 25062 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1445 0 : else { return String(id + " must be an integer\n"); }
1446 : }
1447 350 : else { return String(""); }
1448 : }
1449 :
1450 : // Read Float from Record
1451 39499 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1452 : {
1453 39499 : if( rec.isDefined( id ) )
1454 : {
1455 36704 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1456 36704 : { rec.get( id , val ); return String(""); }
1457 0 : else { return String(id + " must be a float\n"); }
1458 : }
1459 2795 : else { return String(""); }
1460 : }
1461 :
1462 : // Read Bool from Record
1463 49415 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1464 : {
1465 49415 : if( rec.isDefined( id ) )
1466 : {
1467 42104 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1468 0 : else { return String(id + " must be a bool\n"); }
1469 : }
1470 7311 : else{ return String(""); }
1471 : }
1472 :
1473 : // Read Vector<Int> from Record
1474 899 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
1475 : {
1476 899 : if( rec.isDefined( id ) )
1477 : {
1478 899 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1479 899 : { 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 4369 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1493 : {
1494 4369 : if( rec.isDefined( id ) )
1495 : {
1496 4369 : if( rec.dataType( id )==TpArrayFloat )
1497 : {
1498 2574 : 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 1795 : else if ( rec.dataType( id ) ==TpArrayDouble )
1512 : {
1513 170 : Vector<Double> vec; rec.get( id, vec );
1514 170 : val.resize(vec.nelements());
1515 510 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1516 170 : return String("");
1517 170 : }
1518 1625 : else if ( rec.dataType( id ) ==TpArrayInt )
1519 : {
1520 49 : Vector<Int> vec; rec.get( id, vec );
1521 49 : val.resize(vec.nelements());
1522 255 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1523 49 : return String("");
1524 49 : }
1525 1576 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1526 : {
1527 1576 : Vector<Bool> vec; rec.get( id, vec );
1528 1576 : 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 1576 : }
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 4109 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1539 : {
1540 4109 : if( rec.isDefined( id ) )
1541 : {
1542 4109 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1543 4106 : { rec.get( id , val ); return String(""); }
1544 3 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1545 : {
1546 3 : Vector<Bool> vec; rec.get( id, vec );
1547 3 : if( vec.nelements()==0 ){val.resize(0); return String("");}
1548 0 : else{ return String(id + " must be a vector of strings.\n"); }
1549 3 : }
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 11818 : 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 11818 : 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 2148 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1570 : {
1571 : try
1572 : {
1573 2148 : String tmpRF, tmpRA, tmpDEC;
1574 2148 : std::istringstream iss(instr);
1575 2148 : iss >> tmpRF >> tmpRA >> tmpDEC;
1576 2148 : 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 2148 : casacore::Quantity tmpQRA;
1583 2148 : casacore::Quantity tmpQDEC;
1584 2148 : casacore::Quantity::read(tmpQRA, tmpRA);
1585 2148 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1586 :
1587 2148 : 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 2148 : Bool status = MDirection::getType(theRF, tmpRF);
1598 2148 : if (!status) {
1599 2 : throw AipsError();
1600 : }
1601 2146 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1602 2146 : return String("");
1603 2158 : }
1604 2 : catch(AipsError &x)
1605 : {
1606 4 : return String("Error in converting '" + instr + "' to MDirection. Need format : 'J2000 19:59:28.500 +40.44.01.50'\n");
1607 2 : }
1608 : }
1609 :
1610 : // Read Quantity from a Record string
1611 9239 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1612 : {
1613 9239 : if( rec.isDefined( id ) )
1614 : {
1615 8291 : if( rec.dataType( id )==TpString )
1616 8291 : { 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 948 : else{ return String(""); }
1620 : }
1621 :
1622 : // Read MDirection from a Record string
1623 3096 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1624 : {
1625 3096 : if( rec.isDefined( id ) )
1626 : {
1627 2148 : if( rec.dataType( id )==TpString )
1628 2148 : { 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 948 : else{ return String(""); }
1632 : }
1633 :
1634 : // Convert MDirection to String
1635 1798 : String SynthesisParams::MDirectionToString(MDirection val) const
1636 : {
1637 1798 : MVDirection mvpc( val.getAngle() );
1638 1798 : MVAngle ra = mvpc.get()(0);
1639 1798 : MVAngle dec = mvpc.get()(1);
1640 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1641 5394 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1642 1798 : }
1643 :
1644 : // Convert Quantity to String
1645 6058 : String SynthesisParams::QuantityToString(Quantity val) const
1646 : {
1647 6058 : 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 6058 : ss.precision(std::numeric_limits<double>::digits10+2);
1652 6058 : ss << val;
1653 18174 : 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 6058 : }
1659 :
1660 : // Convert Record contains Quantity or Measure quantities to String
1661 49 : String SynthesisParams::recordQMToString(const Record &rec) const
1662 : {
1663 : Double val;
1664 49 : String unit;
1665 49 : if ( rec.isDefined("m0") )
1666 : {
1667 28 : Record subrec = rec.subRecord("m0");
1668 28 : subrec.get("value",val);
1669 28 : subrec.get("unit",unit);
1670 28 : }
1671 21 : else if (rec.isDefined("value") )
1672 : {
1673 21 : rec.get("value",val);
1674 21 : rec.get("unit",unit);
1675 : }
1676 147 : return String::toString(val) + unit;
1677 49 : }
1678 :
1679 :
1680 : /////////////////////// Selection Parameters
1681 :
1682 5112 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1683 : {
1684 5112 : setDefaults();
1685 5112 : }
1686 :
1687 5122 : SynthesisParamsSelect::~SynthesisParamsSelect()
1688 : {
1689 5122 : }
1690 :
1691 11 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1692 11 : operator=(other);
1693 11 : }
1694 :
1695 2037 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1696 2037 : if(this!=&other) {
1697 2037 : msname=other.msname;
1698 2037 : spw=other.spw;
1699 2037 : freqbeg=other.freqbeg;
1700 2037 : freqend=other.freqend;
1701 2037 : freqframe=other.freqframe;
1702 2037 : field=other.field;
1703 2037 : antenna=other.antenna;
1704 2037 : timestr=other.timestr;
1705 2037 : scan=other.scan;
1706 2037 : obs=other.obs;
1707 2037 : state=other.state;
1708 2037 : uvdist=other.uvdist;
1709 2037 : taql=other.taql;
1710 2037 : usescratch=other.usescratch;
1711 2037 : readonly=other.readonly;
1712 2037 : incrmodel=other.incrmodel;
1713 2037 : datacolumn=other.datacolumn;
1714 :
1715 : }
1716 2037 : return *this;
1717 : }
1718 :
1719 3086 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1720 : {
1721 3086 : setDefaults();
1722 :
1723 3086 : String err("");
1724 :
1725 : try
1726 : {
1727 :
1728 3086 : err += readVal( inrec, String("msname"), msname );
1729 :
1730 3086 : err += readVal( inrec, String("readonly"), readonly );
1731 3086 : err += readVal( inrec, String("usescratch"), usescratch );
1732 :
1733 : // override with entries from savemodel.
1734 3086 : String savemodel("");
1735 3086 : err += readVal( inrec, String("savemodel"), savemodel );
1736 3086 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1737 1961 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1738 1944 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1739 :
1740 3086 : err += readVal( inrec, String("incrmodel"), incrmodel );
1741 :
1742 3086 : err += readVal( inrec, String("spw"), spw );
1743 :
1744 : /// -------------------------------------------------------------------------------------
1745 : // Why are these params here ? Repeats of defineimage.
1746 3086 : err += readVal( inrec, String("freqbeg"), freqbeg);
1747 3086 : err += readVal( inrec, String("freqend"), freqend);
1748 :
1749 3086 : String freqframestr( MFrequency::showType(freqframe) );
1750 3086 : err += readVal( inrec, String("outframe"), freqframestr);
1751 3086 : if( ! MFrequency::getType(freqframe, freqframestr) )
1752 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1753 : /// -------------------------------------------------------------------------------------
1754 :
1755 3086 : err += readVal( inrec, String("field"),field);
1756 3086 : err += readVal( inrec, String("antenna"),antenna);
1757 3086 : err += readVal( inrec, String("timestr"),timestr);
1758 3086 : err += readVal( inrec, String("scan"),scan);
1759 3086 : err += readVal( inrec, String("obs"),obs);
1760 3086 : err += readVal( inrec, String("state"),state);
1761 3086 : err += readVal( inrec, String("uvdist"),uvdist);
1762 3086 : err += readVal( inrec, String("taql"),taql);
1763 :
1764 3086 : err += readVal( inrec, String("datacolumn"),datacolumn);
1765 :
1766 3086 : err += verify();
1767 :
1768 3086 : }
1769 0 : catch(AipsError &x)
1770 : {
1771 0 : err = err + x.getMesg() + "\n";
1772 0 : }
1773 :
1774 3086 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1775 :
1776 3086 : }
1777 :
1778 3086 : String SynthesisParamsSelect::verify() const
1779 : {
1780 3086 : String err;
1781 : // Does the MS exist on disk.
1782 3086 : Directory thems( msname );
1783 3086 : if( thems.exists() )
1784 : {
1785 : // Is it readable ?
1786 3086 : if( ! thems.isReadable() )
1787 0 : { err += "MS " + msname + " is not readable.\n"; }
1788 : // Depending on 'readonly', is the MS writable ?
1789 3086 : 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 6172 : return err;
1796 3086 : }
1797 :
1798 8198 : void SynthesisParamsSelect::setDefaults()
1799 : {
1800 8198 : msname="";
1801 8198 : spw="";
1802 8198 : freqbeg="";
1803 8198 : freqend="";
1804 8198 : MFrequency::getType(freqframe,"LSRK");
1805 8198 : field="";
1806 8198 : antenna="";
1807 8198 : timestr="";
1808 8198 : scan="";
1809 8198 : obs="";
1810 8198 : state="";
1811 8198 : uvdist="";
1812 8198 : taql="";
1813 8198 : usescratch=false;
1814 8198 : readonly=true;
1815 8198 : incrmodel=false;
1816 8198 : datacolumn="corrected";
1817 8198 : }
1818 :
1819 2086 : Record SynthesisParamsSelect::toRecord()const
1820 : {
1821 2086 : Record selpar;
1822 2086 : selpar.define("msname",msname);
1823 2086 : selpar.define("spw",spw);
1824 2086 : selpar.define("freqbeg",freqbeg);
1825 2086 : selpar.define("freqend",freqend);
1826 2086 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1827 : //looks like fromRecord looks for outframe !
1828 2086 : selpar.define("outframe", MFrequency::showType(freqframe));
1829 2086 : selpar.define("field",field);
1830 2086 : selpar.define("antenna",antenna);
1831 2086 : selpar.define("timestr",timestr);
1832 2086 : selpar.define("scan",scan);
1833 2086 : selpar.define("obs",obs);
1834 2086 : selpar.define("state",state);
1835 2086 : selpar.define("uvdist",uvdist);
1836 2086 : selpar.define("taql",taql);
1837 2086 : selpar.define("usescratch",usescratch);
1838 2086 : selpar.define("readonly",readonly);
1839 2086 : selpar.define("incrmodel",incrmodel);
1840 2086 : selpar.define("datacolumn",datacolumn);
1841 :
1842 2086 : return selpar;
1843 0 : }
1844 :
1845 :
1846 : /////////////////////// Image Parameters
1847 :
1848 5551 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1849 : {
1850 5551 : setDefaults();
1851 5551 : }
1852 :
1853 5550 : SynthesisParamsImage::~SynthesisParamsImage()
1854 : {
1855 5550 : }
1856 :
1857 3719 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1858 3719 : if(this != &other){
1859 3719 : imageName=other.imageName;
1860 3719 : stokes=other.stokes;
1861 3719 : startModel.resize(); startModel=other.startModel;
1862 3719 : imsize.resize(); imsize=other.imsize;
1863 3719 : cellsize.resize(); cellsize=other.cellsize;
1864 3719 : projection=other.projection;
1865 3719 : useNCP=other.useNCP;
1866 3719 : phaseCenter=other.phaseCenter;
1867 3719 : phaseCenterFieldId=other.phaseCenterFieldId;
1868 3719 : obslocation=other.obslocation;
1869 3719 : pseudoi=other.pseudoi;
1870 3719 : nchan=other.nchan;
1871 3719 : nTaylorTerms=other.nTaylorTerms;
1872 3719 : chanStart=other.chanStart;
1873 3719 : chanStep=other.chanStep;
1874 3719 : freqStart=other.freqStart;
1875 3719 : freqStep=other.freqStep;
1876 3719 : refFreq=other.refFreq;
1877 3719 : velStart=other.velStart;
1878 3719 : velStep=other.velStep;
1879 3719 : freqFrame=other.freqFrame;
1880 3719 : mFreqStart=other.mFreqStart;
1881 3719 : mFreqStep=other.mFreqStep;
1882 3719 : mVelStart=other.mVelStart;
1883 3719 : mVelStep=other.mVelStep;
1884 3719 : restFreq.resize(); restFreq=other.restFreq;
1885 3719 : start=other.start;
1886 3719 : step=other.step;
1887 3719 : frame=other.frame;
1888 3719 : veltype=other.veltype;
1889 3719 : mode=other.mode;
1890 3719 : reffreq=other.reffreq;
1891 3719 : sysvel=other.sysvel;
1892 3719 : sysvelframe=other.sysvelframe;
1893 3719 : sysvelvalue=other.sysvelvalue;
1894 3719 : qmframe=other.qmframe;
1895 3719 : mveltype=other.mveltype;
1896 3719 : tststr=other.tststr;
1897 3719 : startRecord=other.startRecord;
1898 3719 : stepRecord=other.stepRecord;
1899 3719 : reffreqRecord=other.reffreqRecord;
1900 3719 : sysvelRecord=other.sysvelRecord;
1901 3719 : restfreqRecord=other.restfreqRecord;
1902 3719 : csysRecord=other.csysRecord;
1903 3719 : csys=other.csys;
1904 3719 : imshape.resize(); imshape=other.imshape;
1905 3719 : freqFrameValid=other.freqFrameValid;
1906 3719 : overwrite=other.overwrite;
1907 3719 : deconvolver=other.deconvolver;
1908 3719 : distance=other.distance;
1909 3719 : trackDir=other.trackDir;
1910 3719 : trackSource=other.trackSource;
1911 3719 : movingSource=other.movingSource;
1912 :
1913 :
1914 :
1915 : }
1916 :
1917 3719 : return *this;
1918 :
1919 : }
1920 :
1921 1687 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec) {
1922 1687 : fromRecord(inrec, False);
1923 1686 : }
1924 :
1925 1847 : void SynthesisParamsImage::fromRecord(const casacore::Record &inrec,
1926 : const casacore::Bool isSingleDish) {
1927 1847 : setDefaults();
1928 1847 : String err("");
1929 3694 : LogIO os(LogOrigin("SynthesisParamsImage", "fromRecord"));
1930 :
1931 : try {
1932 1847 : err += readVal(inrec, String("imagename"), imageName);
1933 :
1934 : //// imsize
1935 1847 : if (inrec.isDefined("imsize")) {
1936 1847 : DataType tp = inrec.dataType("imsize");
1937 1847 : if (tp == TpInt) { // A single integer for both dimensions.
1938 : Int npix;
1939 742 : inrec.get("imsize", npix);
1940 742 : imsize.resize(2);
1941 742 : imsize.set(npix);
1942 1105 : } else if (tp == TpArrayInt) { // An integer array : [ nx ] or [ nx, ny ]
1943 1105 : Vector<Int> ims;
1944 1105 : inrec.get("imsize", ims);
1945 1105 : if (ims.nelements() == 1) { // [ nx ]
1946 47 : imsize.set(ims[0]);
1947 1058 : } else if (ims.nelements() == 2) { // [ nx, ny ]
1948 1058 : imsize[0] = ims[0];
1949 1058 : 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 1105 : } 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 1847 : if (inrec.isDefined("cell")) {
1960 : try {
1961 1847 : DataType tp = inrec.dataType("cell");
1962 1847 : if (tp == TpInt ||
1963 1847 : tp == TpFloat ||
1964 : tp == TpDouble) {
1965 11 : Double cell = inrec.asDouble("cell");
1966 11 : cellsize.set(Quantity(cell, "arcsec"));
1967 1847 : } else if (tp == TpArrayInt ||
1968 1836 : tp == TpArrayFloat ||
1969 : tp == TpArrayDouble) {
1970 1 : Vector<Double> cells;
1971 1 : inrec.get("cell", cells);
1972 1 : if (cells.nelements() == 1) { // [ cellx ]
1973 0 : cellsize.set(Quantity(cells[0], "arcsec"));
1974 1 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1975 1 : cellsize[0] = Quantity(cells[0], "arcsec");
1976 1 : 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 1836 : } else if (tp == TpString) {
1981 835 : String cell;
1982 835 : inrec.get("cell", cell);
1983 835 : Quantity qcell;
1984 835 : err += stringToQuantity(cell, qcell);
1985 835 : cellsize.set(qcell);
1986 1835 : } else if (tp == TpArrayString) {
1987 1000 : Array<String> cells;
1988 1000 : inrec.get("cell", cells);
1989 1000 : Vector<String> vcells(cells);
1990 1000 : if (cells.nelements() == 1) { // [ cellx ]
1991 12 : Quantity qcell;
1992 12 : err += stringToQuantity(vcells[0], qcell);
1993 12 : cellsize.set(qcell);
1994 1000 : } else if (cells.nelements() == 2) { // [ cellx, celly ]
1995 988 : err += stringToQuantity(vcells[0], cellsize[0]);
1996 988 : 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 1000 : } 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 1847 : err += readVal(inrec, String("stokes"), stokes);
2010 1847 : if (stokes.matches("pseudoI")) {
2011 1 : stokes = "I";
2012 1 : pseudoi = true;
2013 : } else {
2014 1846 : pseudoi = false;
2015 : }
2016 :
2017 : /// PseudoI
2018 :
2019 : ////nchan
2020 1847 : err += readVal(inrec, String("nchan"), nchan);
2021 :
2022 : /// phaseCenter (as a string) . // Add INT support later.
2023 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2024 1847 : if (inrec.isDefined("phasecenter")) {
2025 1847 : String pcent("");
2026 1847 : if (inrec.dataType("phasecenter") == TpString) {
2027 1821 : inrec.get("phasecenter", pcent);
2028 1821 : if (pcent.length() > 0) { // if it's zero length, it means 'figure out from first field in MS'.
2029 1249 : err += readVal(inrec, String("phasecenter"), phaseCenter);
2030 1249 : 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 1249 : err += readVal(inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2034 : } else {
2035 572 : phaseCenterFieldId = -2; // deal with this later in buildCoordinateSystem to assign the first selected field
2036 : }
2037 : }
2038 1847 : if (inrec.dataType("phasecenter") == TpInt ||
2039 5515 : inrec.dataType("phasecenter") == TpFloat ||
2040 3668 : inrec.dataType("phasecenter") == TpDouble) {
2041 : // This will override the previous setting to 0 if the phaseCenter string is zero length.
2042 26 : err += readVal(inrec, String("phasecenter"), phaseCenterFieldId);
2043 : }
2044 1873 : if ((inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter") != TpInt
2045 1873 : && inrec.dataType("phasecenter") != TpFloat && inrec.dataType("phasecenter") != TpDouble)) {
2046 0 : err += String("Cannot set phasecenter. Please specify a string or int\n");
2047 : }
2048 1847 : }
2049 :
2050 : //// Projection
2051 1847 : if (inrec.isDefined("projection")) {
2052 1847 : if (inrec.dataType("projection") == TpString) {
2053 1847 : String pstr;
2054 1847 : inrec.get("projection", pstr);
2055 : try {
2056 1847 : if (pstr.matches("NCP")) {
2057 1 : pstr = "SIN";
2058 1 : useNCP = true;
2059 : }
2060 1847 : projection = Projection::type(pstr);
2061 0 : } catch (AipsError &x) {
2062 0 : err += String("Invalid projection code : " + pstr);
2063 0 : }
2064 1847 : } else {
2065 0 : err += "projection must be a string\n";
2066 : }
2067 : } //projection
2068 :
2069 : // Frequency frame stuff.
2070 1847 : err += readVal(inrec, String("specmode"), mode);
2071 : // Alias for 'mfs' is 'cont'
2072 1847 : if (mode == "cont") mode = "mfs";
2073 :
2074 1847 : err += readVal(inrec, String("outframe"), frame);
2075 1847 : qmframe = "";
2076 : // mveltype is only set when start/step is given in mdoppler
2077 1847 : mveltype = "";
2078 : //start
2079 1847 : String startType("");
2080 1847 : String widthType("");
2081 1847 : if (inrec.isDefined("start")) {
2082 1847 : if (inrec.dataType("start") == TpInt) {
2083 239 : err += readVal(inrec, String("start"), chanStart);
2084 239 : start = String::toString(chanStart);
2085 239 : startType = "chan";
2086 1608 : } else if (inrec.dataType("start") == TpString) {
2087 1573 : err += readVal(inrec, String("start"), start);
2088 1573 : if (start.contains("Hz")) {
2089 157 : stringToQuantity(start, freqStart);
2090 157 : startType = "freq";
2091 1416 : } else if (start.contains("m/s")) {
2092 46 : stringToQuantity(start, velStart);
2093 46 : startType = "vel";
2094 : }
2095 35 : } 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 35 : startRecord = inrec.subRecord("start");
2100 35 : if (startRecord.isDefined("m0")) {
2101 : //must be a measure
2102 21 : String mtype;
2103 21 : startRecord.get("type", mtype);
2104 21 : if (mtype == "frequency") {
2105 : //mfrequency
2106 7 : startRecord.get(String("refer"), qmframe);
2107 7 : 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 7 : start = recordQMToString(startRecord);
2112 7 : stringToQuantity(start, freqStart);
2113 7 : startType = "freq";
2114 14 : } else if (mtype == "radialvelocity") {
2115 : //mradialvelocity
2116 7 : startRecord.get(String("refer"), qmframe);
2117 7 : if (frame != "" && frame != qmframe) {
2118 : // should emit warning to the logger
2119 7 : cerr << "The frame in start:" << qmframe << " Override frame=" << frame << endl;
2120 : }
2121 7 : start = recordQMToString(startRecord);
2122 7 : stringToQuantity(start, velStart);
2123 7 : startType = "vel";
2124 7 : } else if (mtype == "doppler") {
2125 : //use veltype in mdoppler
2126 : //start = MDopToVelString(startRecord);
2127 7 : start = recordQMToString(startRecord);
2128 7 : stringToQuantity(start, velStart);
2129 7 : startRecord.get(String("refer"), mveltype);
2130 7 : mveltype.downcase();
2131 7 : startType = "vel";
2132 : }
2133 21 : } else {
2134 14 : start = recordQMToString(startRecord);
2135 14 : if (start.contains("Hz")) {
2136 7 : stringToQuantity(start, freqStart);
2137 7 : startType = "freq";
2138 7 : } else if (start.contains("m/s")) {
2139 7 : stringToQuantity(start, velStart);
2140 7 : 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 1847 : if (inrec.isDefined("width")) {
2152 1847 : if (inrec.dataType("width") == TpInt) {
2153 206 : err += readVal(inrec, String("width"), chanStep);
2154 206 : step = String::toString(chanStep);
2155 206 : widthType = "chan";
2156 1641 : } else if (inrec.dataType("width") == TpString) {
2157 1627 : err += readVal(inrec, String("width"), step);
2158 1627 : if (step.contains("Hz")) {
2159 142 : stringToQuantity(step, freqStep);
2160 142 : widthType = "freq";
2161 1485 : } else if (step.contains("m/s")) {
2162 50 : stringToQuantity(step, velStep);
2163 50 : widthType = "vel";
2164 : }
2165 14 : } 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 14 : stepRecord = inrec.subRecord("width");
2170 14 : if (stepRecord.isDefined("m0")) {
2171 : //must be a measure
2172 7 : String mtype;
2173 7 : stepRecord.get("type", mtype);
2174 7 : 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 7 : } else if (mtype == "radialvelocity") {
2185 : //mradialvelocity
2186 7 : stepRecord.get(String("refer"), qmframe);
2187 7 : 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 7 : step = recordQMToString(stepRecord);
2192 7 : stringToQuantity(step, velStep);
2193 7 : 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 7 : } else {
2203 7 : step = recordQMToString(stepRecord);
2204 7 : if (step.contains("Hz")) {
2205 0 : stringToQuantity(step, freqStep);
2206 0 : widthType = "freq";
2207 7 : } else if (step.contains("m/s")) {
2208 7 : stringToQuantity(step, velStep);
2209 7 : 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 1847 : 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 1847 : if (inrec.isDefined("reffreq")) {
2226 1847 : if (inrec.dataType("reffreq") == TpString) {
2227 1847 : 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 1847 : err += readVal(inrec, String("veltype"), veltype);
2254 1847 : veltype = mveltype != "" ? mveltype : veltype;
2255 : // sysvel (String, Quantity)
2256 1847 : if (inrec.isDefined("sysvel")) {
2257 1847 : if (inrec.dataType("sysvel") == TpString) {
2258 1847 : 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 1847 : err += readVal(inrec, String("sysvelframe"), sysvelframe);
2270 :
2271 : // rest frequencies (record or vector of Strings)
2272 1847 : if (inrec.isDefined("restfreq")) {
2273 1847 : Vector<String> rfreqs(0);
2274 1847 : Record restfreqSubRecord;
2275 1847 : 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 1847 : } else if (inrec.dataType("restfreq") == TpArrayString) {
2286 967 : err += readVal(inrec, String("restfreq"), rfreqs);
2287 880 : } else if (inrec.dataType("restfreq") == TpString) {
2288 7 : rfreqs.resize(1);
2289 7 : err += readVal(inrec, String("restfreq"), rfreqs[0]);
2290 : }
2291 1847 : restFreq.resize(rfreqs.nelements());
2292 2107 : for (uInt fr = 0; fr < rfreqs.nelements(); fr++) {
2293 260 : err += stringToQuantity(rfreqs[fr], restFreq[fr]);
2294 : }
2295 1847 : } // if def restfreq
2296 :
2297 : // optional - coordsys, imshape
2298 : // if exist use them. May need a consistency check with the rest of impars?
2299 1847 : if (inrec.isDefined("csys")) {
2300 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
2301 899 : if (inrec.dataType("csys") == TpRecord) {
2302 : //csysRecord = inrec.subRecord("csys");
2303 899 : csysRecord.defineRecord("coordsys", inrec.subRecord("csys"));
2304 : }
2305 899 : if (inrec.isDefined("imshape")) {
2306 899 : if (inrec.dataType("imshape") == TpArrayInt) {
2307 899 : 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 1847 : String freqframestr = (frame != "" && qmframe != "") ? qmframe : frame;
2318 1847 : if (frame != "" && !MFrequency::getType(freqFrame, freqframestr)) {
2319 1 : err += "Invalid Frequency Frame " + freqframestr;
2320 : }
2321 1847 : err += readVal(inrec, String("restart"), overwrite);
2322 :
2323 1847 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2324 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2325 1847 : if (inrec.isDefined("startmodel")) {
2326 1847 : if (inrec.dataType("startmodel") == TpString) {
2327 941 : String onemodel;
2328 941 : err += readVal(inrec, String("startmodel"), onemodel);
2329 941 : if (onemodel.length() > 0) {
2330 16 : startModel.resize(1);
2331 16 : startModel[0] = onemodel;
2332 : } else {
2333 925 : startModel.resize();
2334 : }
2335 2754 : } else if (inrec.dataType("startmodel") == TpArrayString ||
2336 907 : inrec.dataType("startmodel") == TpArrayBool) {
2337 906 : 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 1847 : err += readVal(inrec, String("nterms"), nTaylorTerms);
2344 1847 : err += readVal(inrec, String("deconvolver"), deconvolver);
2345 :
2346 : // Force nchan=1 for anything other than cube modes...
2347 1847 : if (mode == "mfs") nchan = 1;
2348 : //read obslocation
2349 1847 : if (inrec.isDefined("obslocation_rec")) {
2350 899 : String errorobs;
2351 899 : const Record obsrec = inrec.asRecord("obslocation_rec");
2352 899 : MeasureHolder mh;
2353 899 : if (!mh.fromRecord(errorobs, obsrec)) {
2354 0 : err += errorobs;
2355 : }
2356 899 : obslocation = mh.asMPosition();
2357 899 : }
2358 1847 : err += verify(isSingleDish);
2359 1847 : }
2360 0 : catch (AipsError &x) {
2361 0 : err = err + x.getMesg() + "\n";
2362 0 : }
2363 :
2364 1847 : if (err.length() > 0) throw(AipsError("Invalid Image Parameter set : " + err));
2365 1851 : }
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 1847 : String SynthesisParamsImage::verify(const casacore::Bool isSingleDish) const
2401 : {
2402 1847 : LogIO os;
2403 1847 : String err;
2404 :
2405 1847 : if(imageName=="") {err += "Please supply an image name\n";}
2406 :
2407 1847 : if(imsize.nelements() != 2){err += "imsize must be a vector of 2 Ints\n";}
2408 1847 : if(cellsize.nelements() != 2) {err += "cellsize must be a vector of 2 Quantities\n";}
2409 1847 : if(cellsize[0].getValue() == 0.0 || cellsize[1].getValue() == 0.0){
2410 1 : 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 1847 : if((mode=="mfs") && nchan > 1){
2418 0 : err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";
2419 : }
2420 :
2421 2073 : if(! stokes.matches("I") && ! stokes.matches("Q") &&
2422 215 : ! stokes.matches("U") && ! stokes.matches("V") &&
2423 153 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2424 153 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2425 151 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2426 140 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2427 138 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2428 2201 : ! stokes.matches("IQU") && ! stokes.matches("IUV") &&
2429 128 : ! 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 1847 : if(startModel.nelements()>0){
2438 22 : if(deconvolver!="mtmfs"){
2439 :
2440 16 : if( startModel.nelements()!=1 ){
2441 0 : err += String("Only one startmodel image is allowed.\n");
2442 : }
2443 : else{
2444 32 : File fp(imageName+String(".model"));
2445 16 : 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 16 : }
2447 : }
2448 : else{// mtmfs
2449 12 : File fp(imageName+String(".model.tt0"));
2450 6 : 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 6 : }
2456 :
2457 : // Check that startmodel exists on disk !
2458 49 : for(uInt ss=0;ss<startModel.nelements();ss++){
2459 27 : File fp( startModel[ss] );
2460 27 : if(!fp.exists()){
2461 0 : err += "Startmodel " + startModel[ss] + " cannot be found on disk.";
2462 : }
2463 27 : }
2464 : }
2465 :
2466 : /// Check imsize for efficiency.
2467 1847 : Int imxnew = SynthesisUtilMethods::getOptimumSize(imsize[0]);
2468 1847 : Int imynew = SynthesisUtilMethods::getOptimumSize(imsize[1]);
2469 1847 : if((imxnew != imsize[0] || imynew != imsize[1]) && isSingleDish == False){
2470 18 : LogIO os(LogOrigin("SynthesisParamsImage","fromRecord",WHERE));
2471 9 : if( imxnew != imsize[0] ) {
2472 2 : os << LogIO::WARN << "imsize with "+String::toString(imsize[0])+
2473 2 : " pixels is not an efficient imagesize. Try "+
2474 3 : String::toString(imxnew)+" instead." << LogIO::POST;
2475 : }
2476 9 : if( imsize[0] != imsize[1] && imynew != imsize[1] ){
2477 16 : os << LogIO::WARN << "imsize with "+String::toString(imsize[1])+
2478 16 : " pixels is not an efficient imagesize. Try "+
2479 24 : String::toString(imynew)+" instead." << LogIO::POST;
2480 : }
2481 9 : }
2482 3694 : return err;
2483 1847 : }// 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 7398 : void SynthesisParamsImage::setDefaults()
2496 : {
2497 : // Image definition parameters
2498 7398 : imageName = String("");
2499 7398 : imsize.resize(2); imsize.set(100);
2500 7398 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2501 7398 : stokes="I";
2502 7398 : phaseCenter=MDirection();
2503 7398 : phaseCenterFieldId=-1;
2504 7398 : projection=Projection::SIN;
2505 7398 : useNCP=false;
2506 7398 : startModel=Vector<String>(0);
2507 7398 : freqFrameValid=True;
2508 7398 : overwrite=false;
2509 : // PseudoI
2510 7398 : pseudoi=false;
2511 :
2512 : // Spectral coordinates
2513 7398 : nchan=1;
2514 7398 : mode="mfs";
2515 7398 : start="";
2516 7398 : step="";
2517 7398 : chanStart=0;
2518 7398 : 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 7398 : freqStart=Quantity(0,"");
2524 7398 : freqStep=Quantity(0,"");
2525 7398 : velStart=Quantity(0,"");
2526 7398 : velStep=Quantity(0,"");
2527 7398 : veltype=String("radio");
2528 7398 : restFreq.resize(0);
2529 7398 : refFreq = Quantity(0,"Hz");
2530 7398 : frame = "";
2531 7398 : freqFrame=MFrequency::LSRK;
2532 7398 : sysvel="";
2533 7398 : sysvelframe="";
2534 7398 : sysvelvalue=Quantity(0.0,"m/s");
2535 7398 : nTaylorTerms=1;
2536 7398 : deconvolver="hogbom";
2537 : ///csysRecord=Record();
2538 7398 : }
2539 :
2540 899 : Record SynthesisParamsImage::toRecord() const
2541 : {
2542 899 : Record impar;
2543 899 : impar.define("imagename", imageName);
2544 899 : impar.define("imsize", imsize);
2545 899 : Vector<String> cells(2);
2546 899 : cells[0] = QuantityToString( cellsize[0] );
2547 899 : cells[1] = QuantityToString( cellsize[1] );
2548 899 : impar.define("cell", cells );
2549 899 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2550 899 : else{impar.define("stokes", stokes);}
2551 899 : impar.define("nchan", nchan);
2552 899 : impar.define("nterms", nTaylorTerms);
2553 899 : impar.define("deconvolver",deconvolver);
2554 899 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2555 899 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2556 899 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2557 :
2558 899 : impar.define("specmode", mode );
2559 : // start and step can be one of these types
2560 899 : if( start!="" )
2561 : {
2562 299 : if( !start.contains("Hz") && !start.contains("m/s") &&
2563 65 : String::toInt(start) == chanStart )
2564 : {
2565 65 : impar.define("start",chanStart);
2566 : }
2567 169 : else if( startRecord.nfields() > 0 )
2568 : {
2569 25 : impar.defineRecord("start", startRecord );
2570 : }
2571 : else
2572 : {
2573 144 : impar.define("start",start);
2574 : }
2575 : }
2576 : else {
2577 665 : impar.define("start", start );
2578 : }
2579 899 : if( step!="" )
2580 : {
2581 228 : if( !step.contains("Hz") && !step.contains("m/s") &&
2582 41 : String::toInt(step) == chanStep )
2583 : {
2584 41 : impar.define("width", chanStep);
2585 : }
2586 146 : else if( stepRecord.nfields() > 0 )
2587 : {
2588 10 : impar.defineRecord("width",stepRecord);
2589 : }
2590 : else
2591 : {
2592 136 : impar.define("width",step);
2593 : }
2594 : }
2595 : else
2596 : {
2597 712 : 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 899 : impar.define("veltype", veltype);
2606 899 : if (restfreqRecord.nfields() != 0 )
2607 : {
2608 0 : impar.defineRecord("restfreq", restfreqRecord);
2609 : }
2610 : else
2611 : {
2612 899 : Vector<String> rfs( restFreq.nelements() );
2613 1084 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2614 899 : impar.define("restfreq", rfs);
2615 899 : }
2616 : //impar.define("reffreq", QuantityToString(refFreq));
2617 : //reffreq
2618 899 : if( reffreqRecord.nfields() != 0 )
2619 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2620 : else
2621 899 : { impar.define("reffreq",reffreq); }
2622 : //impar.define("reffreq", reffreq );
2623 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2624 899 : impar.define("outframe", frame );
2625 : //sysvel
2626 899 : if( sysvelRecord.nfields() != 0 )
2627 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2628 : else
2629 899 : { impar.define("sysvel", sysvel );}
2630 899 : impar.define("sysvelframe", sysvelframe );
2631 :
2632 899 : impar.define("restart",overwrite );
2633 899 : impar.define("freqframevalid", freqFrameValid);
2634 899 : impar.define("startmodel", startModel );
2635 :
2636 899 : if( csysRecord.isDefined("coordsys") )
2637 : {
2638 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2639 899 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2640 899 : impar.define("imshape", imshape);
2641 : }
2642 : // else cout << " NO CSYS INFO to write to output record " << endl;
2643 : ///Now save obslocation
2644 899 : Record tmprec;
2645 899 : String err;
2646 899 : MeasureHolder mh(obslocation);
2647 899 : if(mh.toRecord(err, tmprec)){
2648 899 : impar.defineRecord("obslocation_rec", tmprec);
2649 : }
2650 : else{
2651 0 : throw(AipsError("failed to save obslocation to record"));
2652 : }
2653 1798 : return impar;
2654 899 : }
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 945 : 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 945 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2674 945 : vi2.originChunks();
2675 945 : vi2.origin();
2676 : /// This version uses the new vi2/vb2
2677 : // get the first ms for multiple MSes
2678 : //MeasurementSet msobj=vi2.ms();
2679 945 : Int fld=vb->fieldId()(0);
2680 :
2681 : //handling first ms only
2682 945 : Double gfreqmax=-1.0;
2683 945 : Double gdatafend=-1.0;
2684 945 : Double gdatafstart=1e14;
2685 945 : Double gfreqmin=1e14;
2686 945 : Vector<Int> spwids0;
2687 945 : Int j=0;
2688 945 : Int minfmsid=0;
2689 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2690 945 : Double imStartFreq=getCubeImageStartFreq();
2691 945 : std::vector<Int> sourceMsWithStartFreq;
2692 :
2693 :
2694 1954 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2695 : //auto forMS0=chansel.find(0);
2696 1009 : map<Int, Vector<Int> > spwsels=forMS0->second;
2697 1009 : Int nspws=spwsels.size();
2698 1009 : Vector<Int> spwids(nspws);
2699 1009 : Vector<Int> nChannels(nspws);
2700 1009 : Vector<Int> firstChannels(nspws);
2701 : //Vector<Int> channelIncrement(nspws);
2702 :
2703 1009 : Int k=0;
2704 2315 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2705 1306 : spwids[k]=it->first;
2706 1306 : nChannels[k]=(it->second)[0];
2707 1306 : firstChannels[k]=(it->second)[1];
2708 : }
2709 1009 : if(j==0) {
2710 945 : spwids0.resize();
2711 945 : 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 1009 : Double freqmin=0, freqmax=0;
2723 1009 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2724 :
2725 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2726 1009 : 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 1009 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2732 : //cerr << "after " << datafstart << " " << datafend << endl;
2733 1009 : if((datafstart > datafend) || !status)
2734 0 : throw(AipsError("spw selection failed"));
2735 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2736 :
2737 1009 : if (mode=="cubedata") {
2738 2 : freqmin = datafstart;
2739 2 : freqmax = datafend;
2740 : }
2741 1007 : else if(mode == "cubesource"){
2742 5 : if(!trackSource){
2743 0 : throw(AipsError("Cannot be in cubesource without tracking a moving source"));
2744 : }
2745 5 : String ephemtab(movingSource);
2746 5 : if(movingSource=="TRACKFIELD"){
2747 3 : Int fieldID=MSColumns(*mss[j]).fieldId()(0);
2748 3 : ephemtab=Path(MSColumns(*mss[j]).field().ephemPath(fieldID)).absoluteName();
2749 : }
2750 5 : MEpoch refep=MSColumns(*mss[j]).timeMeas()(0);
2751 5 : Quantity refsysvel;
2752 5 : MSUtil::getFreqRangeAndRefFreqShift(freqmin,freqmax,refsysvel, refep, spwids,firstChannels, nChannels, *mss[j], ephemtab, trackDir, true);
2753 5 : if(j==0)
2754 5 : 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 5 : }
2762 : else {
2763 : //VisBufferUtil::getFreqRange(freqmin,freqmax, vi2, freqFrameValid? freqFrame:MFrequency::REST );
2764 : //cerr << "before " << freqmin << " " << freqmax << endl;
2765 2004 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2766 2004 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2767 : //cerr << "after " << freqmin << " " << freqmax << endl;
2768 : }
2769 :
2770 :
2771 :
2772 :
2773 1009 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2774 1009 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2775 1009 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2776 1009 : 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 1009 : if(imStartFreq > 0.0 && imStartFreq >= freqmin && imStartFreq <= freqmax){
2783 72 : if(mode != "cubesource"){
2784 72 : minfmsid=j;
2785 72 : spwids0.resize();
2786 72 : spwids0=spwids;
2787 72 : vi2.originChunks();
2788 72 : vi2.origin();
2789 74 : while(vb->msId() != j && vi2.moreChunks() ){
2790 2 : vi2.nextChunk();
2791 2 : vi2.origin();
2792 : }
2793 72 : fld=vb->fieldId()(0);
2794 :
2795 : }
2796 : else{
2797 0 : sourceMsWithStartFreq.push_back(j);
2798 : }
2799 : }
2800 :
2801 1009 : }
2802 945 : 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 945 : MeasurementSet msobj = *mss[minfmsid];
2809 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2810 1890 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2811 947 : }
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 945 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2844 : MeasurementSet& msobj,
2845 : Vector<Int> spwids, Int fld,
2846 : Double freqmin, Double freqmax,
2847 : Double datafstart, Double datafend )
2848 : {
2849 1890 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2850 :
2851 945 : CoordinateSystem csys;
2852 945 : if( csysRecord.nfields()!=0 )
2853 : {
2854 : //use cysRecord
2855 1 : Record subRec1;
2856 : // cout<<"USE THE EXISTING CSYS +++++++++++++++++"<<endl;
2857 1 : CoordinateSystem *csysptr = CoordinateSystem::restore(csysRecord,"coordsys");
2858 : //csys = *csysptr;
2859 : //CoordinateSystem csys(*csysptr);
2860 1 : csys = *csysptr;
2861 :
2862 1 : }
2863 : else {
2864 944 : MSColumns msc(msobj);
2865 944 : String telescop = msc.observation().telescopeName()(0);
2866 944 : MEpoch obsEpoch = msc.timeMeas()(0);
2867 944 : MPosition obsPosition;
2868 944 : 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 944 : MDirection phaseCenterToUse = phaseCenter;
2880 :
2881 944 : if( phaseCenterFieldId != -1 )
2882 : {
2883 598 : MSFieldColumns msfield(msobj.field());
2884 598 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2885 : {
2886 572 : if(trackSource){
2887 14 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2888 : }
2889 : else{
2890 558 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2891 : }
2892 : }
2893 : else
2894 : {
2895 26 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2896 : }
2897 598 : }
2898 : // Setup Phase center if it is specified only by field id.
2899 :
2900 : /////////////////// Direction Coordinates
2901 944 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2902 : // Normalize correctly
2903 944 : MVAngle ra=mvPhaseCenter.get()(0);
2904 944 : ra(0.0);
2905 :
2906 944 : MVAngle dec=mvPhaseCenter.get()(1);
2907 944 : Vector<Double> refCoord(2);
2908 944 : refCoord(0)=ra.get().getValue();
2909 944 : refCoord(1)=dec;
2910 944 : Vector<Double> refPixel(2);
2911 944 : refPixel(0) = Double(imsize[0]/2);
2912 944 : refPixel(1) = Double(imsize[1]/2);
2913 :
2914 944 : Vector<Double> deltas(2);
2915 944 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2916 944 : deltas(1)=cellsize[1].get("rad").getValue();
2917 944 : Matrix<Double> xform(2,2);
2918 944 : xform=0.0;xform.diagonal()=1.0;
2919 :
2920 944 : Vector<Double> projparams(2);
2921 944 : projparams = 0.0;
2922 944 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2923 944 : Projection projTo( projection.type(), projparams );
2924 :
2925 : DirectionCoordinate
2926 944 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2927 : // projection,
2928 : projTo,
2929 2832 : refCoord(0), refCoord(1),
2930 2832 : deltas(0), deltas(1),
2931 : xform,
2932 1888 : 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 944 : obslocation=obsPosition;
2940 944 : ObsInfo myobsinfo;
2941 944 : myobsinfo.setTelescope(telescop);
2942 944 : myobsinfo.setPointingCenter(mvPhaseCenter);
2943 944 : myobsinfo.setObsDate(obsEpoch);
2944 944 : 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 944 : if(spwids.nelements()==0)
2960 : {
2961 0 : Int nspw=msc.spectralWindow().nrow();
2962 0 : spwids.resize(nspw);
2963 0 : indgen(spwids);
2964 : }
2965 944 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2966 944 : Vector<Double> dataChanFreq, dataChanWidth;
2967 944 : std::vector<std::vector<Int> > averageWhichChan;
2968 944 : std::vector<std::vector<Int> > averageWhichSPW;
2969 944 : std::vector<std::vector<Double> > averageChanFrac;
2970 :
2971 944 : if(spwids.nelements()==1)
2972 : {
2973 787 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2974 787 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
2975 : }
2976 : else
2977 : {
2978 157 : 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 944 : Double minDataFreq = min(dataChanFreq);
2985 944 : if(start=="" && minDataFreq < datafstart ) {
2986 : // limit data chan freq vector for default start case with channel selection
2987 : Int chanStart, chanEnd;
2988 6 : Int lochan = 0;
2989 6 : Int nDataChan = dataChanFreq.nelements();
2990 6 : Int hichan = nDataChan-1;
2991 : Double diff_fmin, diff_fmax;
2992 6 : Bool ascending = dataChanFreq[nDataChan-1] - dataChanFreq[0] > 0;
2993 90 : for(Int ichan = 0; ichan < nDataChan; ichan++)
2994 : {
2995 84 : diff_fmin = dataChanFreq[ichan] - datafstart;
2996 84 : diff_fmax = datafend - dataChanFreq[ichan];
2997 : // freqmin and freqmax should corresponds to the channel edges
2998 84 : if(ascending)
2999 : {
3000 :
3001 84 : if( diff_fmin > 0 && diff_fmin <= dataChanWidth[ichan]/2. )
3002 : {
3003 6 : lochan = ichan;
3004 : }
3005 78 : else if(diff_fmax > 0 && diff_fmax <= dataChanWidth[ichan]/2. )
3006 : {
3007 4 : 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 6 : chanStart = lochan;
3023 6 : chanEnd = hichan;
3024 6 : if (lochan > hichan)
3025 : {
3026 0 : chanStart=hichan;
3027 0 : chanEnd=lochan;
3028 : }
3029 6 : Vector<Double> tempChanFreq = dataChanFreq(Slice(chanStart,chanEnd-chanStart+1,1));
3030 6 : Vector<Double> tempChanWidth = dataChanWidth(Slice(chanStart,chanEnd-chanStart+1,1));
3031 6 : dataChanFreq.resize(tempChanFreq.nelements());
3032 6 : dataChanWidth.resize(tempChanWidth.nelements());
3033 6 : dataChanFreq = tempChanFreq;
3034 6 : dataChanWidth = tempChanWidth;
3035 6 : }
3036 944 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3037 944 : String cubemode;
3038 944 : if ( qrestfreq.getValue("Hz")==0 )
3039 : {
3040 870 : MSDopplerUtil msdoppler(msobj);
3041 870 : Vector<Double> restfreqvec;
3042 870 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3043 870 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3044 870 : if ( qrestfreq.getValue("Hz")==0 and mode!="mfs" )
3045 : {
3046 90 : cubemode = findSpecMode(mode);
3047 90 : 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 90 : Double provisional_restfreq = (datafend+datafstart)/2.0;
3052 90 : 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 90 : << 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 870 : }
3062 : Double refPix;
3063 944 : Vector<Double> chanFreq;
3064 944 : Vector<Double> chanFreqStep;
3065 944 : String specmode;
3066 :
3067 944 : if(mode=="cubesource"){
3068 5 : MDoppler mdop(sysvelvalue, MDoppler::RELATIVISTIC);
3069 5 : dataChanFreq=mdop.shiftFrequency(dataChanFreq);
3070 5 : dataChanWidth=mdop.shiftFrequency(dataChanWidth);
3071 5 : 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 5 : }
3076 :
3077 944 : 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 943 : Bool nonLinearFreq(false);
3083 943 : String veltype_p=veltype;
3084 943 : veltype_p.upcase();
3085 2825 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3086 2825 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3087 : {
3088 2 : nonLinearFreq= true;
3089 : }
3090 :
3091 943 : SpectralCoordinate mySpectral;
3092 : Double stepf;
3093 943 : if(!nonLinearFreq)
3094 : //TODO: velocity mode default start case (use last channels?)
3095 : {
3096 941 : Double startf=chanFreq[0];
3097 : //Double stepf=chanFreqStep[0];
3098 941 : if(chanFreq.nelements()==1)
3099 : {
3100 586 : stepf=chanFreqStep[0];
3101 : }
3102 : else
3103 : {
3104 355 : stepf=chanFreq[1]-chanFreq[0];
3105 : }
3106 941 : Double restf=qrestfreq.getValue("Hz");
3107 : //stepf=9e8;
3108 941 : 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 941 : if(mode=="cubedata")
3112 : {
3113 : // mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3114 4 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST?
3115 : MFrequency::REST : MFrequency::Undefined,
3116 2 : startf, stepf, refPix, restf);
3117 : }
3118 939 : 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 10 : mySpectral = SpectralCoordinate(MFrequency::REST,
3124 5 : startf, stepf, refPix, restf);
3125 : }
3126 : else
3127 : {
3128 1868 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3129 934 : startf, stepf, refPix, restf);
3130 : }
3131 : }
3132 : else
3133 : { // nonlinear freq coords - use tabular setting
3134 : // once NOFRAME is implemented do this
3135 2 : 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 2 : else if (mode=="cubesource")
3143 : {
3144 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
3145 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3146 : }
3147 : else
3148 : {
3149 4 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3150 6 : 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 943 : uInt nrestfreq = restFreq.nelements();
3159 943 : 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 943 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3177 943 : if(whichStokes.nelements()==0)
3178 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3179 943 : StokesCoordinate myStokes(whichStokes);
3180 :
3181 : ////////////////// Build Full coordinate system.
3182 :
3183 : //CoordinateSystem csys;
3184 943 : csys.addCoordinate(myRaDec);
3185 943 : csys.addCoordinate(myStokes);
3186 943 : csys.addCoordinate(mySpectral);
3187 943 : csys.setObsInfo(myobsinfo);
3188 :
3189 : //store back csys to impars record
3190 : //cerr<<"save csys to csysRecord..."<<endl;
3191 943 : if(csysRecord.isDefined("coordsys"))
3192 0 : csysRecord.removeField("coordsys");
3193 943 : csys.save(csysRecord,"coordsys");
3194 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3195 : // imshape
3196 943 : imshape.resize(4);
3197 943 : imshape[0] = imsize[0];
3198 943 : imshape[1] = imsize[0];
3199 943 : imshape[2] = whichStokes.nelements();
3200 943 : 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 969 : } // end of else when coordsys record is not defined...
3234 :
3235 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3236 :
3237 1888 : return csys;
3238 946 : }
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 14 : MDirection SynthesisParamsImage::getMovingSourceDir(const MeasurementSet& ms, const MEpoch& refEp, const MPosition& obsposition, const MDirection::Types outframe){
3490 14 : MDirection outdir;
3491 14 : String ephemtab(movingSource);
3492 14 : if(movingSource=="TRACKFIELD"){
3493 6 : Int fieldID=MSColumns(ms).fieldId()(0);
3494 6 : ephemtab=Path(MSColumns(ms).field().ephemPath(fieldID)).absoluteName();
3495 : }
3496 14 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
3497 14 : 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 14 : MeasFrame mframe(refEp, obsposition);
3500 14 : MDirection::Ref outref1(MDirection::AZEL, mframe);
3501 14 : MDirection::Ref outref(outframe, mframe);
3502 14 : 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 14 : if (Table::isReadable(ephemtab)){
3508 12 : MeasComet mcomet(Path(ephemtab).absoluteName());
3509 12 : mframe.set(mcomet);
3510 12 : tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
3511 12 : }
3512 2 : else if (planetType >=MDirection::MERCURY && planetType <MDirection::COMET){
3513 2 : tmpazel=MDirection::Convert(trackDir, outref1)();
3514 : }
3515 14 : outdir=MDirection::Convert(tmpazel, outref)();
3516 :
3517 28 : return outdir;
3518 14 : }
3519 :
3520 944 : 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 944 : String inStart, inStep;
3531 944 : specmode = findSpecMode(mode);
3532 944 : String freqframe;
3533 944 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3534 1888 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3535 :
3536 944 : refPix=0.0;
3537 944 : Bool descendingfreq(false);
3538 944 : Bool descendingoutfreq(false);
3539 :
3540 944 : if( mode.contains("cube") )
3541 : {
3542 479 : String restfreq=QuantityToString(qrestfreq);
3543 : // use frame from input start or width in MFreaquency or MRadialVelocity
3544 479 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3545 : // emit warning here if qmframe is used
3546 : //
3547 479 : inStart = start;
3548 479 : inStep = step;
3549 479 : if( specmode=="channel" )
3550 : {
3551 398 : inStart = String::toString(chanStart);
3552 398 : inStep = String::toString(chanStep);
3553 : // negative step -> descending channel indices
3554 398 : if (inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3555 : // input frame is the data frame
3556 : //freqframe = MFrequency::showType(dataFrame);
3557 : }
3558 81 : 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 49 : if(inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3569 : }
3570 32 : 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 32 : if (!inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3582 : }
3583 :
3584 479 : if (inStep=='0') inStep="";
3585 :
3586 479 : MRadialVelocity mSysVel;
3587 479 : Quantity qVel;
3588 : MRadialVelocity::Types mRef;
3589 479 : if(mode!="cubesource")
3590 : {
3591 :
3592 :
3593 474 : 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 5 : freqframe=MFrequency::showType(dataFrame);
3603 5 : 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 5 : { mSysVel=MRadialVelocity();}
3610 : }
3611 : // cubedata mode: input start, step are those of the input data frame
3612 479 : if ( mode=="cubedata" )
3613 : {
3614 2 : freqframe=MFrequency::showType(dataFrame);
3615 2 : 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 479 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3626 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3627 958 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3628 479 : << LogIO::POST;
3629 479 : ostringstream ostr;
3630 479 : ostr << " phaseCenter='" << phaseCenter;
3631 479 : os << String(ostr)<<"' ";
3632 :
3633 : Double dummy; // dummy variable - weightScale is not used here
3634 958 : 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 479 : veltype,
3651 : verbose,
3652 : mSysVel
3653 : );
3654 :
3655 479 : if( nchan==-1 )
3656 : {
3657 280 : nchan = chanFreq.nelements();
3658 280 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3659 : }
3660 :
3661 479 : if (!rst) {
3662 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
3663 1 : << LogIO::EXCEPTION;
3664 0 : return false;
3665 : }
3666 : os << LogIO::DEBUG1
3667 478 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3668 956 : << LogIO::POST;
3669 :
3670 478 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3671 14 : descendingoutfreq = true;
3672 : }
3673 :
3674 : //if (descendingfreq && !descendingoutfreq) {
3675 876 : if ((specmode=="channel" && descendingfreq==1)
3676 876 : || (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 31 : std::vector<Double> stlchanfreq;
3682 31 : chanFreq.tovector(stlchanfreq);
3683 31 : std::reverse(stlchanfreq.begin(),stlchanfreq.end());
3684 31 : chanFreq=Vector<Double>(stlchanfreq);
3685 31 : chanFreqStep=-chanFreqStep;
3686 31 : }
3687 482 : }
3688 465 : else if ( mode=="mfs" ) {
3689 465 : chanFreq.resize(1);
3690 465 : chanFreqStep.resize(1);
3691 : //chanFreqStep[0] = freqmax - freqmin;
3692 465 : Double freqmean = (freqmin + freqmax)/2;
3693 465 : if (refFreq.getValue("Hz")==0) {
3694 409 : chanFreq[0] = freqmean;
3695 409 : refPix = 0.0;
3696 409 : chanFreqStep[0] = freqmax - freqmin;
3697 : }
3698 : else {
3699 56 : chanFreq[0] = refFreq.getValue("Hz");
3700 : // Set the new reffreq to be the refPix (CAS-9518)
3701 56 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3702 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3703 56 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3704 : }
3705 :
3706 465 : if( nchan==-1 ) nchan=1;
3707 465 : if( qrestfreq.getValue("Hz")==0.0 ) {
3708 48 : restFreq.resize(1);
3709 48 : 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 943 : return true;
3719 :
3720 947 : }//getImFreq
3721 : /////////////////////////
3722 945 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3723 945 : Double inStartFreq=-1.0;
3724 945 : String checkspecmode("");
3725 945 : if(mode.contains("cube")) {
3726 480 : checkspecmode = findSpecMode(mode);
3727 : }
3728 945 : if(checkspecmode!="") {
3729 480 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3730 480 : if(checkspecmode=="channel") {
3731 399 : inStartFreq=-1.0;
3732 : }
3733 : else {
3734 81 : if(checkspecmode=="frequency") {
3735 49 : inStartFreq = freqStart.get("Hz").getValue();
3736 : }
3737 32 : else if(checkspecmode=="velocity") {
3738 : MDoppler::Types DopType;
3739 32 : MDoppler::getType(DopType, veltype);
3740 32 : MDoppler mdop(velStart,DopType);
3741 32 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3742 32 : inStartFreq = MFrequency::fromDoppler(mdop, qrestfreq.getValue(Unit("Hz")), mfreqframe).getValue();
3743 32 : }
3744 : }
3745 : }
3746 :
3747 945 : return inStartFreq;
3748 :
3749 945 : }
3750 :
3751 1514 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3752 : {
3753 1514 : String specmode;
3754 1514 : specmode="channel";
3755 1514 : 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 2054 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3761 1005 : !( velStep.getValue()==0.0)) {
3762 64 : specmode="velocity";
3763 : }
3764 1875 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3765 890 : !(freqStep.getValue()==0.0)) {
3766 103 : specmode="frequency";
3767 : }
3768 : }
3769 1514 : return specmode;
3770 0 : }
3771 :
3772 :
3773 2830 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3774 : {
3775 2830 : Vector<Int> whichStokes(0);
3776 3301 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3777 567 : stokes=="RR" ||stokes=="LL" ||
3778 3301 : stokes=="XX" || stokes=="YY" ) {
3779 2647 : whichStokes.resize(1);
3780 2647 : whichStokes(0)=Stokes::type(stokes);
3781 : }
3782 519 : else if(stokes=="IV" || stokes=="IQ" ||
3783 498 : stokes=="RRLL" || stokes=="XXYY" ||
3784 513 : stokes=="QU" || stokes=="UV"){
3785 33 : whichStokes.resize(2);
3786 :
3787 33 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3788 18 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3789 18 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3790 18 : else if(stokes=="XXYY"){whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
3791 12 : 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 150 : 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 150 : else if(stokes=="IQUV"){
3804 150 : whichStokes.resize(4);
3805 150 : whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;
3806 150 : whichStokes(2)=Stokes::U; whichStokes(3)=Stokes::V;
3807 : }
3808 :
3809 2830 : return whichStokes;
3810 0 : }// decidenpolplanes
3811 :
3812 1887 : IPosition SynthesisParamsImage::shp() const
3813 : {
3814 1887 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3815 :
3816 1887 : if( imsize[0]<=0 || imsize[1]<=0 || nStokes<=0 || nchan<=0 )
3817 : {
3818 1 : 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 3772 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3822 : }
3823 :
3824 943 : Record SynthesisParamsImage::getcsys() const
3825 : {
3826 943 : 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 5551 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3883 : {
3884 5551 : setDefaults();
3885 5551 : }
3886 :
3887 5550 : SynthesisParamsGrid::~SynthesisParamsGrid()
3888 : {
3889 5550 : }
3890 :
3891 :
3892 1847 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3893 : {
3894 1847 : setDefaults();
3895 :
3896 1847 : String err("");
3897 :
3898 : try {
3899 1847 : err += readVal( inrec, String("imagename"), imageName );
3900 :
3901 : // FTMachine parameters
3902 1847 : err += readVal( inrec, String("gridder"), gridder );
3903 1847 : err += readVal( inrec, String("padding"), padding );
3904 1847 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3905 1847 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3906 1847 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3907 1847 : err += readVal( inrec, String("convfunc"), convFunc );
3908 :
3909 1847 : err += readVal( inrec, String("vptable"), vpTable );
3910 :
3911 :
3912 :
3913 : // convert 'gridder' to 'ftmachine' and 'mtype'
3914 1847 : ftmachine = "gridft";
3915 1847 : mType = "default";
3916 1847 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
3917 1270 : ftmachine = "gridft";
3918 : }
3919 623 : else if ( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) &&
3920 46 : (wprojplanes>1 || wprojplanes==-1) ) {
3921 45 : ftmachine = "wprojectft";
3922 : }
3923 : //facetting alone use gridft
3924 532 : else if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes==1))
3925 1 : {ftmachine=="gridft";}
3926 :
3927 531 : else if (gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" ) {
3928 222 : ftmachine = "mosaicft";
3929 : }
3930 :
3931 309 : else if (gridder=="imagemosaic") {
3932 0 : mType = "imagemosaic";
3933 0 : if (wprojplanes>1 || wprojplanes==-1) {
3934 0 : ftmachine = "wprojectft";
3935 : }
3936 : }
3937 :
3938 309 : else if (gridder=="awproject" || gridder=="awprojectft" || gridder=="awp") {
3939 93 : ftmachine = "awprojectft";
3940 : }
3941 :
3942 216 : else if (gridder=="singledish") {
3943 :
3944 160 : ftmachine = "sd";
3945 : }
3946 : else{
3947 56 : ftmachine=gridder;
3948 56 : ftmachine.downcase();
3949 :
3950 : }
3951 :
3952 1847 : String deconvolver;
3953 1847 : err += readVal( inrec, String("deconvolver"), deconvolver );
3954 1847 : if (deconvolver=="mtmfs") {
3955 122 : mType = "multiterm"; // Takes precedence over imagemosaic
3956 : }
3957 :
3958 : // facets
3959 1847 : err += readVal( inrec, String("facets"), facets );
3960 : // chanchunks
3961 1847 : err += readVal( inrec, String("chanchunks"), chanchunks );
3962 :
3963 : // Spectral interpolation
3964 1847 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
3965 : // Track moving source ?
3966 1847 : err += readVal( inrec, String("distance"), distance );
3967 1847 : err += readVal( inrec, String("tracksource"), trackSource );
3968 1847 : err += readVal( inrec, String("trackdir"), trackDir );
3969 :
3970 : // The extra params for WB-AWP
3971 1847 : err += readVal( inrec, String("aterm"), aTermOn );
3972 1847 : err += readVal( inrec, String("psterm"), psTermOn );
3973 1847 : err += readVal( inrec, String("mterm"), mTermOn );
3974 1847 : err += readVal( inrec, String("wbawp"), wbAWP );
3975 1847 : err += readVal( inrec, String("cfcache"), cfCache );
3976 1847 : err += readVal( inrec, String("usepointing"), usePointing );
3977 1847 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3978 1847 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3979 1847 : err += readVal( inrec, String("conjbeams"), conjBeams );
3980 1847 : err += readVal( inrec, String("computepastep"), computePAStep );
3981 1847 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3982 :
3983 : // The extra params for single-dish
3984 1847 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3985 1847 : err += readVal( inrec, String("convertfirst"), convertFirst );
3986 1847 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3987 1847 : err += readVal( inrec, String("convsupport"), convSupport );
3988 1847 : err += readVal( inrec, String("truncate"), truncateSize );
3989 1847 : err += readVal( inrec, String("gwidth"), gwidth );
3990 1847 : err += readVal( inrec, String("jwidth"), jwidth );
3991 1847 : err += readVal( inrec, String("minweight"), minWeight );
3992 1847 : 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 1847 : if (ftmachine=="awprojectft" && cfCache=="") {
3998 0 : cfCache = imageName + ".cf";
3999 : }
4000 :
4001 1847 : if ( ftmachine=="awprojectft" &&
4002 1881 : usePointing==True &&
4003 34 : pointingOffsetSigDev.nelements() != 2 ) {
4004 : // Set the default to a large value so that it behaves like CASA 5.6's usepointing=True.
4005 8 : pointingOffsetSigDev.resize(2);
4006 8 : pointingOffsetSigDev[0] = 600.0;
4007 8 : pointingOffsetSigDev[1] = 600.0;
4008 : }
4009 :
4010 1847 : err += verify();
4011 :
4012 1847 : } catch(AipsError &x) {
4013 0 : err = err + x.getMesg() + "\n";
4014 0 : }
4015 :
4016 1847 : if (err.length()>0) {
4017 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4018 : }
4019 :
4020 1847 : }
4021 :
4022 1847 : String SynthesisParamsGrid::verify() const
4023 : {
4024 1847 : 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 1847 : if ( imageName == "" ) {
4032 0 : err += "Please supply an image name\n";
4033 : }
4034 :
4035 1107 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4036 2847 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4037 2583 : (ftmachine != "mawprojectft") &&
4038 160 : (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 3787 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4049 1940 : ( 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 1847 : 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 1847 : 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 1847 : 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 1847 : 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 1847 : 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 1847 : 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 1881 : if ( ftmachine == "awprojectft" and usePointing == True and
4107 34 : 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 1847 : if ( ftmachine == "sd" ) {
4115 320 : if ( convertFirst != "always" and
4116 320 : 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 1847 : return err;
4124 0 : }
4125 :
4126 7398 : void SynthesisParamsGrid::setDefaults()
4127 : {
4128 7398 : imageName="";
4129 : // FTMachine parameters
4130 : //ftmachine="GridFT";
4131 7398 : ftmachine="gridft";
4132 7398 : gridder=ftmachine;
4133 7398 : padding=1.2;
4134 7398 : useAutoCorr=false;
4135 7398 : useDoublePrec=true;
4136 7398 : wprojplanes=1;
4137 7398 : convFunc="SF";
4138 7398 : vpTable="";
4139 :
4140 : // facets
4141 7398 : facets=1;
4142 :
4143 : // chanchunks
4144 7398 : chanchunks=1;
4145 :
4146 : // Spectral Axis interpolation
4147 7398 : interpolation=String("nearest");
4148 :
4149 : //mosaic use pointing
4150 7398 : usePointing=false;
4151 : // Moving phase center ?
4152 7398 : distance=Quantity(0,"m");
4153 7398 : trackSource=false;
4154 7398 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4155 :
4156 : // The extra params for WB-AWP
4157 7398 : aTermOn = true;
4158 7398 : psTermOn = true;
4159 7398 : mTermOn = false;
4160 7398 : wbAWP = true;
4161 7398 : cfCache = "";
4162 7398 : usePointing = false;
4163 7398 : pointingOffsetSigDev.resize(0);
4164 : // pointingOffsetSigDev.set(30.0);
4165 7398 : doPBCorr = true;
4166 7398 : conjBeams = true;
4167 7398 : computePAStep=360.0;
4168 7398 : rotatePAStep=5.0;
4169 :
4170 : // extra params for single-dish
4171 7398 : pointingDirCol = "";
4172 7398 : convertFirst = "never";
4173 7398 : skyPosThreshold = 0.0;
4174 7398 : convSupport = -1;
4175 7398 : truncateSize = Quantity(-1.0);
4176 7398 : gwidth = Quantity(-1.0);
4177 7398 : jwidth = Quantity(-1.0);
4178 7398 : minWeight = 0.0;
4179 7398 : clipMinMax = False;
4180 :
4181 : // Mapper type
4182 7398 : mType = String("default");
4183 :
4184 7398 : }
4185 :
4186 899 : Record SynthesisParamsGrid::toRecord() const
4187 : {
4188 899 : Record gridpar;
4189 :
4190 899 : gridpar.define("imagename", imageName);
4191 : // FTMachine params
4192 899 : gridpar.define("padding", padding);
4193 899 : gridpar.define("useautocorr",useAutoCorr );
4194 899 : gridpar.define("usedoubleprec", useDoublePrec);
4195 899 : gridpar.define("wprojplanes", wprojplanes);
4196 899 : gridpar.define("convfunc", convFunc);
4197 899 : gridpar.define("vptable", vpTable);
4198 :
4199 899 : gridpar.define("facets", facets);
4200 899 : gridpar.define("chanchunks", chanchunks);
4201 :
4202 899 : gridpar.define("interpolation",interpolation);
4203 :
4204 899 : gridpar.define("distance", QuantityToString(distance));
4205 899 : gridpar.define("tracksource", trackSource);
4206 899 : gridpar.define("trackdir", MDirectionToString( trackDir ));
4207 :
4208 899 : gridpar.define("aterm",aTermOn );
4209 899 : gridpar.define("psterm",psTermOn );
4210 899 : gridpar.define("mterm",mTermOn );
4211 899 : gridpar.define("wbawp", wbAWP);
4212 899 : gridpar.define("cfcache", cfCache);
4213 899 : gridpar.define("usepointing",usePointing );
4214 899 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
4215 899 : gridpar.define("dopbcorr", doPBCorr);
4216 899 : gridpar.define("conjbeams",conjBeams );
4217 899 : gridpar.define("computepastep", computePAStep);
4218 899 : gridpar.define("rotatepastep", rotatePAStep);
4219 :
4220 899 : gridpar.define("pointingcolumntouse", pointingDirCol );
4221 899 : gridpar.define("convertfirst", convertFirst );
4222 899 : gridpar.define("skyposthreshold", skyPosThreshold );
4223 899 : gridpar.define("convsupport", convSupport );
4224 899 : gridpar.define("truncate", QuantityToString(truncateSize) );
4225 899 : gridpar.define("gwidth", QuantityToString(gwidth) );
4226 899 : gridpar.define("jwidth", QuantityToString(jwidth) );
4227 899 : gridpar.define("minweight", minWeight );
4228 899 : gridpar.define("clipminmax", clipMinMax );
4229 :
4230 899 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
4231 : /// else gridpar.define("deconvolver","singleterm");
4232 :
4233 899 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
4234 899 : else gridpar.define("gridder", gridder);
4235 :
4236 : // gridpar.define("mtype", mType);
4237 :
4238 899 : return gridpar;
4239 0 : }
4240 :
4241 :
4242 :
4243 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
4244 :
4245 : /////////////////////// Deconvolver Parameters
4246 :
4247 2526 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4248 : {
4249 2526 : setDefaults();
4250 2526 : }
4251 :
4252 2526 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4253 : {
4254 2526 : }
4255 :
4256 :
4257 2522 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4258 : {
4259 2522 : setDefaults();
4260 :
4261 2522 : String err("");
4262 :
4263 : try
4264 : {
4265 :
4266 2522 : err += readVal( inrec, String("imagename"), imageName );
4267 2522 : 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 2522 : if( inrec.isDefined("startmodel") )
4273 : {
4274 2522 : if( inrec.dataType("startmodel")==TpString )
4275 : {
4276 838 : String onemodel;
4277 838 : err += readVal( inrec, String("startmodel"), onemodel );
4278 838 : if( onemodel.length()>0 )
4279 : {
4280 22 : startModel.resize(1);
4281 22 : startModel[0] = onemodel;
4282 : }
4283 816 : else {startModel.resize();}
4284 838 : }
4285 3370 : else if( inrec.dataType("startmodel")==TpArrayString ||
4286 1686 : inrec.dataType("startmodel")==TpArrayBool)
4287 : {
4288 1684 : 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 2522 : err += readVal( inrec, String("id"), deconvolverId );
4297 2522 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4298 :
4299 2522 : err += readVal( inrec, String("scales"), scales );
4300 2522 : err += readVal( inrec, String("scalebias"), scalebias );
4301 :
4302 2522 : err += readVal( inrec, String("usemask"), maskType );
4303 2522 : if( maskType=="auto-thresh" )
4304 : {
4305 0 : autoMaskAlgorithm = "thresh";
4306 : }
4307 2522 : else if( maskType=="auto-thesh2" )
4308 : {
4309 0 : autoMaskAlgorithm = "thresh2";
4310 : }
4311 2522 : else if( maskType=="auto-onebox" )
4312 : {
4313 0 : autoMaskAlgorithm = "onebox";
4314 : }
4315 2522 : else if( maskType=="user" || maskType=="pb" )
4316 : {
4317 2415 : autoMaskAlgorithm = "";
4318 : }
4319 :
4320 :
4321 2522 : if( inrec.isDefined("mask") )
4322 : {
4323 2522 : if( inrec.dataType("mask")==TpString )
4324 : {
4325 1968 : err+= readVal( inrec, String("mask"), maskString );
4326 : }
4327 554 : else if( inrec.dataType("mask")==TpArrayString )
4328 : {
4329 552 : err+= readVal( inrec, String("mask"), maskList );
4330 : }
4331 : }
4332 :
4333 2522 : if( inrec.isDefined("pbmask") )
4334 : {
4335 2522 : err += readVal( inrec, String("pbmask"), pbMask );
4336 : }
4337 2522 : if( inrec.isDefined("maskthreshold") )
4338 : {
4339 2522 : if( inrec.dataType("maskthreshold")==TpString )
4340 : {
4341 2522 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4342 : //deal with the case a string is a float value without unit
4343 2522 : Quantity testThresholdString;
4344 2522 : Quantity::read(testThresholdString,maskThreshold);
4345 2522 : if( testThresholdString.getUnit()=="" )
4346 : {
4347 2522 : if(testThresholdString.getValue()<1.0)
4348 : {
4349 2522 : fracOfPeak = testThresholdString.getValue();
4350 2522 : maskThreshold=String("");
4351 : }
4352 : }
4353 2522 : }
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 2522 : if( inrec.isDefined("maskresolution") )
4372 : {
4373 2522 : if( inrec.dataType("maskresolution")==TpString )
4374 : {
4375 2522 : err += readVal(inrec, String("maskresolution"), maskResolution );
4376 : //deal with the case a string is a float value without unit
4377 2522 : Quantity testResolutionString;
4378 2522 : Quantity::read(testResolutionString,maskResolution);
4379 2522 : if( testResolutionString.getUnit()=="" )
4380 : {
4381 2522 : maskResByBeam = testResolutionString.getValue();
4382 2522 : maskResolution=String("");
4383 : }
4384 2522 : }
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 2522 : if( inrec.isDefined("nmask") )
4398 : {
4399 2522 : if( inrec.dataType("nmask")==TpInt )
4400 : {
4401 2522 : err+= readVal(inrec, String("nmask"), nMask );
4402 : }
4403 : else
4404 : {
4405 0 : err+= "nmask must be an integer\n";
4406 : }
4407 : }
4408 2522 : if( inrec.isDefined("autoadjust") )
4409 : {
4410 1675 : if( inrec.dataType("autoadjust")==TpBool )
4411 : {
4412 1675 : 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 2522 : if (inrec.isDefined("fusedthreshold"))
4421 : {
4422 2522 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4423 2522 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4424 : else
4425 0 : err += "fusedthreshold must be a float or double";
4426 : }
4427 2522 : if (inrec.isDefined("specmode"))
4428 : {
4429 2522 : if(inrec.dataType("specmode") == TpString)
4430 2522 : 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 2522 : if (inrec.isDefined("largestscale"))
4436 : {
4437 2522 : if (inrec.dataType("largestscale") == TpInt)
4438 2522 : 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 2522 : if( inrec.isDefined("sidelobethreshold"))
4444 : {
4445 2522 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4446 : {
4447 2522 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4448 : }
4449 : else
4450 : {
4451 0 : err+= "sidelobethreshold must be a float or double";
4452 : }
4453 : }
4454 :
4455 2522 : if( inrec.isDefined("noisethreshold"))
4456 : {
4457 2522 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4458 : {
4459 2522 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4460 : }
4461 : else
4462 : {
4463 0 : err+= "noisethreshold must be a float or double";
4464 : }
4465 : }
4466 2522 : if( inrec.isDefined("lownoisethreshold"))
4467 : {
4468 2522 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4469 : {
4470 2522 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4471 : }
4472 : else
4473 : {
4474 0 : err+= "lownoisethreshold must be a float or double";
4475 : }
4476 : }
4477 2522 : if( inrec.isDefined("negativethreshold"))
4478 : {
4479 2522 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4480 : {
4481 2522 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4482 : }
4483 : else
4484 : {
4485 0 : err+= "negativethreshold must be a float or double";
4486 : }
4487 : }
4488 2522 : if( inrec.isDefined("smoothfactor"))
4489 : {
4490 2522 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4491 : {
4492 2522 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4493 : }
4494 : else
4495 : {
4496 0 : err+= "smoothfactor must be a float or double";
4497 : }
4498 : }
4499 2522 : if( inrec.isDefined("minbeamfrac"))
4500 : {
4501 2522 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4502 : {
4503 2522 : 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 2522 : if( inrec.isDefined("cutthreshold"))
4517 : {
4518 2522 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4519 : {
4520 2522 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4521 : }
4522 : else {
4523 0 : err+= "cutthreshold must be a float or double";
4524 : }
4525 : }
4526 2522 : if( inrec.isDefined("growiterations"))
4527 : {
4528 2522 : if (inrec.dataType("growiterations")==TpInt) {
4529 2522 : err+= readVal(inrec, String("growiterations"), growIterations );
4530 : }
4531 : else {
4532 0 : err+= "growiterations must be an integer\n";
4533 : }
4534 : }
4535 2522 : if( inrec.isDefined("dogrowprune"))
4536 : {
4537 2522 : if (inrec.dataType("dogrowprune")==TpBool) {
4538 2522 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4539 : }
4540 : else {
4541 0 : err+= "dogrowprune must be a bool\n";
4542 : }
4543 : }
4544 2522 : if( inrec.isDefined("minpercentchange"))
4545 : {
4546 2522 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4547 2522 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4548 : }
4549 : else {
4550 0 : err+= "minpercentchange must be a float or double";
4551 : }
4552 : }
4553 2522 : if( inrec.isDefined("verbose"))
4554 : {
4555 2522 : if (inrec.dataType("verbose")==TpBool ) {
4556 2522 : err+= readVal(inrec, String("verbose"), verbose);
4557 : }
4558 : else {
4559 0 : err+= "verbose must be a bool";
4560 : }
4561 : }
4562 2522 : if( inrec.isDefined("fastnoise"))
4563 : {
4564 2522 : if (inrec.dataType("fastnoise")==TpBool ) {
4565 2522 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4566 : }
4567 : else {
4568 0 : err+= "fastnoise must be a bool";
4569 : }
4570 : }
4571 2522 : if( inrec.isDefined("nsigma") )
4572 : {
4573 2522 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4574 2522 : 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 2522 : if( inrec.isDefined("noRequireSumwt") )
4587 : {
4588 1861 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4589 1861 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4590 : }
4591 : else {
4592 0 : err+= "noRequireSumwt must be a bool";
4593 : }
4594 : }
4595 2522 : if( inrec.isDefined("fullsummary") )
4596 : {
4597 2522 : if (inrec.dataType("fullsummary")==TpBool) {
4598 2522 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4599 : }
4600 : else {
4601 0 : err+= "fullsummary must be a bool";
4602 : }
4603 : }
4604 2522 : if( inrec.isDefined("restoringbeam") )
4605 : {
4606 847 : String errinfo("");
4607 : try {
4608 :
4609 847 : if( inrec.dataType("restoringbeam")==TpString )
4610 : {
4611 16 : 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 16 : if( (! usebeam.matches("common")) && usebeam.length()!=0 )
4616 : {
4617 4 : Quantity bsize;
4618 4 : err += readVal( inrec, String("restoringbeam"), bsize );
4619 4 : restoringbeam.setMajorMinor( bsize, bsize );
4620 4 : usebeam = String("");
4621 4 : }
4622 16 : errinfo = usebeam;
4623 : }
4624 831 : 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 847 : }// if isdefined(restoringbeam)
4658 :
4659 2522 : if( inrec.isDefined("interactive") )
4660 : {
4661 2522 : if( inrec.dataType("interactive")==TpBool )
4662 2522 : {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 2522 : err += verify();
4670 :
4671 : }
4672 0 : catch(AipsError &x)
4673 : {
4674 0 : err = err + x.getMesg() + "\n";
4675 0 : }
4676 :
4677 2522 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4678 :
4679 2522 : }
4680 :
4681 2522 : String SynthesisParamsDeconv::verify() const
4682 : {
4683 2522 : String err;
4684 :
4685 2522 : 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 2522 : if( maskString.length()>0 )
4689 : {
4690 170 : File fp( imageName+".mask" );
4691 170 : 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 170 : }
4693 :
4694 2522 : if( pbMask >= 1.0)
4695 0 : {err += "pbmask must be < 1.0 \n"; }
4696 2522 : else if( pbMask < 0.0)
4697 0 : {err += "pbmask must be a positive value \n"; }
4698 :
4699 2522 : 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 2522 : if ( fracOfPeak >= 1.0)
4708 0 : {err += "fracofpeak must be < 1.0 \n"; }
4709 2522 : else if ( fracOfPeak < 0.0)
4710 0 : {err += "fracofpeak must be a positive value \n"; }
4711 :
4712 2522 : return err;
4713 0 : }
4714 :
4715 5048 : void SynthesisParamsDeconv::setDefaults()
4716 : {
4717 5048 : imageName="";
4718 5048 : algorithm="hogbom";
4719 5048 : startModel=Vector<String>(0);
4720 5048 : deconvolverId=0;
4721 5048 : nTaylorTerms=1;
4722 5048 : scales.resize(1); scales[0]=0.0;
4723 5048 : scalebias=0.6;
4724 5048 : maskType="none";
4725 5048 : maskString="";
4726 5048 : maskList.resize(1); maskList[0]="";
4727 5048 : pbMask=0.0;
4728 5048 : autoMaskAlgorithm="thresh";
4729 5048 : maskThreshold="";
4730 5048 : maskResolution="";
4731 5048 : fracOfPeak=0.0;
4732 5048 : nMask=0;
4733 5048 : interactive=false;
4734 5048 : autoAdjust=False;
4735 5048 : fusedThreshold = 0.0;
4736 5048 : specmode="mfs";
4737 5048 : largestscale = -1;
4738 5048 : }
4739 :
4740 1675 : Record SynthesisParamsDeconv::toRecord() const
4741 : {
4742 1675 : Record decpar;
4743 :
4744 1675 : decpar.define("imagename", imageName);
4745 1675 : decpar.define("deconvolver", algorithm);
4746 1675 : decpar.define("startmodel",startModel);
4747 1675 : decpar.define("id",deconvolverId);
4748 1675 : decpar.define("nterms",nTaylorTerms);
4749 1675 : decpar.define("scales",scales);
4750 1675 : decpar.define("scalebias",scalebias);
4751 1675 : decpar.define("usemask",maskType);
4752 1675 : decpar.define("fusedthreshold", fusedThreshold);
4753 1675 : decpar.define("specmode", specmode);
4754 1675 : decpar.define("largestscale", largestscale);
4755 1675 : if( maskList.nelements()==1 && maskList[0]=="")
4756 : {
4757 1123 : decpar.define("mask",maskString);
4758 : }
4759 : else {
4760 552 : decpar.define("mask",maskList);
4761 : }
4762 1675 : decpar.define("pbmask",pbMask);
4763 1675 : if (fracOfPeak > 0.0)
4764 : {
4765 0 : decpar.define("maskthreshold",fracOfPeak);
4766 : }
4767 : else
4768 : {
4769 1675 : decpar.define("maskthreshold",maskThreshold);
4770 : }
4771 1675 : decpar.define("maskresolution",maskResolution);
4772 1675 : decpar.define("nmask",nMask);
4773 1675 : decpar.define("autoadjust",autoAdjust);
4774 1675 : decpar.define("sidelobethreshold",sidelobeThreshold);
4775 1675 : decpar.define("noisethreshold",noiseThreshold);
4776 1675 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4777 1675 : decpar.define("negativethreshold",negativeThreshold);
4778 1675 : decpar.define("smoothfactor",smoothFactor);
4779 1675 : decpar.define("minbeamfrac",minBeamFrac);
4780 1675 : decpar.define("cutthreshold",cutThreshold);
4781 1675 : decpar.define("growiterations",growIterations);
4782 1675 : decpar.define("dogrowprune",doGrowPrune);
4783 1675 : decpar.define("minpercentchange",minPercentChange);
4784 1675 : decpar.define("verbose", verbose);
4785 1675 : decpar.define("fastnoise", fastnoise);
4786 1675 : decpar.define("interactive",interactive);
4787 1675 : decpar.define("nsigma",nsigma);
4788 1675 : decpar.define("noRequireSumwt",noRequireSumwt);
4789 1675 : decpar.define("fullsummary",fullsummary);
4790 :
4791 1675 : return decpar;
4792 0 : }
4793 :
4794 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4795 :
4796 :
4797 : } //# NAMESPACE CASA - END
4798 :
|
