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