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 0 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 0 : }
89 :
90 0 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 0 : }
93 :
94 0 : Int SynthesisUtilMethods::validate(const VisBuffer& vb)
95 : {
96 0 : Int N=vb.nRow(),M=-1;
97 0 : for(Int i=0;i<N;i++)
98 : {
99 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
100 0 : {M++;break;}
101 : }
102 0 : return M;
103 : }
104 :
105 0 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 0 : Int N=vb.nRows(),M=-1;
108 0 : for(Int i=0;i<N;i++)
109 : {
110 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 0 : {M++;break;}
112 : }
113 0 : return M;
114 : }
115 : // Get the next largest even composite of 2,3,5.
116 : // This is to ensure a 'good' image size for FFTW.
117 : // Translated from gcwrap/scripts/cleanhelper.py : getOptimumSize
118 0 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 0 : Int n=npix;
121 :
122 0 : if( n%2 !=0 ){ n+= 1; }
123 :
124 0 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 0 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 0 : if( fac[k]>5 )
129 : {
130 0 : val = fac[k];
131 0 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 0 : fac[k] = val;
133 : }
134 : }
135 0 : newlarge=product(fac);
136 0 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 0 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 0 : return newlarge;
141 0 : }
142 :
143 : // Return the prime factors of the given number
144 0 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 0 : Vector<uInt> factors;
147 :
148 0 : Int lastresult = n;
149 0 : Int sqlast = int(sqrt(n))+1;
150 :
151 0 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 0 : Int c=2;
153 : while(1)
154 : {
155 0 : if( lastresult == 1 || c > sqlast ) { break; }
156 0 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 0 : if( c>sqlast){ c=lastresult; break; }
160 0 : if( lastresult % c == 0 ) { break; }
161 0 : c += 1;
162 : }
163 0 : factors.resize( factors.nelements()+1, true );
164 0 : factors[ factors.nelements()-1 ] = c;
165 0 : lastresult /= c;
166 : }
167 0 : if( factors.nelements()==0 ) { factors.resize(1);factors[0]=n; }
168 :
169 : //if( douniq ) { factors = unique(factors); }
170 :
171 : /*
172 : /// The Sort::unique isn't working as called below. Results appear to be the
173 : /// same as with the cleanhelper python code, so leaving as is for not. CAS-7889
174 : if( douniq )
175 : {
176 : cout << "Test unique fn on : " << factors << endl;
177 : Sort srt;
178 : Vector<uInt> unvec=factors; uInt nrec;
179 : srt.unique(unvec,nrec);
180 : cout << " Unique : " << unvec << " nr : " << nrec << endl;
181 : }
182 : */
183 :
184 0 : return factors;
185 0 : }
186 :
187 :
188 0 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 0 : if (psfcutoff >=1.0 || psfcutoff<=0.0)
193 : {
194 0 : os << "psfcutoff must be >0 and <1" << LogIO::WARN;
195 0 : return false;
196 : }
197 :
198 0 : std::shared_ptr<SIImageStore> imstore;
199 0 : if( nterms>1 )
200 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 0 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 0 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 0 : imstore->printBeamSet();
210 :
211 0 : imstore->releaseLocks();
212 :
213 0 : return true;
214 0 : }
215 :
216 :
217 :
218 : // Evaluate a polynomial from Taylor coefficients and expand to a Cube.
219 : // Since this funcion is primarily for use with mtmfs_via_cube, this step applies only to
220 : // the multiterm.model.ttx images and the cube.model image cube.
221 0 : Bool SynthesisUtilMethods::taylorCoeffsToCube(const String& cubename,const String& mtname, const Int nterms, const String& reffreq)
222 : {
223 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "taylorCoeffsToCube"));
224 :
225 : // Set up imstores
226 0 : CountedPtr<SIImageStore> cube_imstore;
227 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
228 :
229 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
230 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
231 :
232 : // Check that .model exists.
233 : try{
234 0 : cube_imstore->model();
235 0 : for(Int i=0;i<nterms;i++)
236 0 : mt_imstore->model(i);
237 : }
238 0 : catch(AipsError &x)
239 : {
240 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\nModel images must exist on disk." ));
241 0 : }
242 :
243 : // Get/check shapes
244 0 : IPosition cube_shp( cube_imstore->model()->shape() );
245 0 : IPosition mt_shp( mt_imstore->model(0)->shape() );
246 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
247 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
248 : }
249 :
250 : // Read reference frequency
251 0 : Quantity reffreq_qa;
252 0 : Quantity::read( reffreq_qa, reffreq );
253 0 : Double refval = reffreq_qa.getValue("Hz");
254 :
255 : //cout << "ref freq : " << refval << endl;
256 :
257 : //Get the frequency list for the cube
258 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
259 0 : Vector<Double> freqlist( cube_shp[3] );
260 :
261 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
262 : {
263 0 : if( csys.type(i) == Coordinate::SPECTRAL )
264 : {
265 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
266 :
267 0 : for(Int ch=0;ch<cube_shp[3];ch++)
268 : {
269 : Double freq;
270 0 : Bool ret = speccoord.toWorld( freq, ch );
271 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
272 0 : freqlist[ch] = freq;
273 :
274 : //cout << "freq " << ch << " is " << freq << endl;
275 : }
276 :
277 0 : }
278 : }
279 :
280 :
281 : // Reset the Cube values to zero.
282 : //cube_imstore->model()->set(0.0);
283 :
284 : //For each pol, do the Taylor-to-Cube calculation.
285 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
286 : {
287 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(nterms);
288 0 : for(Int i=0;i<nterms;i++)
289 : {
290 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
291 0 : 0, cube_shp[3],
292 0 : pol, cube_shp[2],
293 0 : *mt_imstore->model(i) );
294 : }
295 :
296 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
297 : {
298 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
299 0 : chan, cube_shp[3],
300 0 : pol, cube_shp[2],
301 0 : *cube_imstore->model() );
302 :
303 0 : Double wt = (freqlist[chan] - refval) / refval;
304 :
305 0 : cube_subim->set(0.0);
306 0 : for(Int tt=0;tt<nterms;tt++)
307 : {
308 0 : Double fac = pow(wt,tt);
309 0 : LatticeExpr<Float> oneterm = LatticeExpr<Float>( *cube_subim + (fac) * (*mt_subims[tt])) ;
310 0 : cube_subim->copyData(oneterm);
311 0 : }
312 :
313 :
314 0 : }//for chan
315 0 : }// for pol
316 :
317 0 : return True;
318 :
319 0 : }//end of func
320 :
321 :
322 : // Calculate the RHS of the Normal equations for a linear least squares fit of a Taylor polynomial (per pixel).
323 : // This function is primarily for use with mtmfs_via_cube, and may be used for the PSF (2nterms-1), the residual (nterms)
324 : // and the primary beam (nterms=1).
325 : // imtype=0 : PSF with 2nterms-1 terms
326 : // imtype=1 : residual with nterms terms
327 : // imtype=2 : pb with 1 term
328 : // imtype=3 : sumwt with 2nterms-1 terms
329 0 : Bool SynthesisUtilMethods::cubeToTaylorSum(const String& cubename,const String& mtname, const Int nterms, const String& reffreq, const Int imtype, const Float pblimit)
330 : {
331 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "cubeToTaylorSum"));
332 :
333 : //cout << "imtype : " << imtype << endl;
334 0 : if(imtype <0 || imtype >3)
335 : {
336 0 : throw( AipsError("cubeToTaylorSum currently only supports 'psf','residual','pb', 'sumwt' options"));
337 : }
338 :
339 : // Set up imstores
340 0 : CountedPtr<SIImageStore> cube_imstore;
341 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
342 :
343 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
344 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, nterms, true, true ));
345 : //For psf avg pb has to be done already
346 : // doing residual too for sensitivity this is independent of beam spectral index removal
347 0 : Float maxPB=1.0;
348 0 : if(imtype < 2){
349 0 : LatticeExprNode elnod( max( *(mt_imstore->pb(0)) ) );
350 0 : maxPB=elnod.getFloat();
351 0 : if(maxPB == 0.0){
352 0 : throw(AipsError("Programmers error: should do tt psf images after making average PB"));
353 :
354 : }
355 :
356 0 : }
357 : // cerr << "imtype " << imtype << " MAX PB " << maxPB << endl;
358 : // If dopsf=True, calculate 2n-1 terms.
359 0 : Int out_nterms=nterms; // for residual
360 0 : if(imtype==0 || imtype==3){out_nterms=2 * nterms - 1;} // the psfs fill the upper triangle of the Hessian with 2 nterms-1 elements. Also sumwt.
361 0 : if(imtype==2){out_nterms=1;} // For the PB, for mtmfs_via_cube, we need only tt0. Later, if we need all terms to calculate PB alpha, then change this to nterms, and add the invHesian math (elsewhere) to later convert the RHS vector into the coefficients.
362 :
363 0 : CountedPtr <ImageInterface<Float> > use_cube, use_mt;
364 : // If dopsf=True, check that .psf cube and mt's exist. If dopsf=False, check residual images.
365 : try{
366 0 : switch(imtype)
367 : {
368 0 : case 0: use_cube=cube_imstore->psf();break;
369 0 : case 1: use_cube=cube_imstore->residual();break;
370 0 : case 2: use_cube=cube_imstore->pb();break;
371 0 : case 3: use_cube=cube_imstore->sumwt();break;
372 : }
373 0 : cube_imstore->sumwt();
374 0 : for(Int i=0;i<out_nterms;i++)
375 : {
376 0 : switch(imtype)
377 : {
378 0 : case 0:mt_imstore->psf(i);break;
379 0 : case 1:mt_imstore->residual(i);break;
380 0 : case 2:mt_imstore->pb(i);break;
381 0 : case 3:mt_imstore->sumwt(i);break;
382 : }
383 : }
384 : }
385 0 : catch(AipsError &x)
386 : {
387 0 : throw( AipsError("Error in reading image : " + x.getMesg() + "\n " + imtype + " images must exist on disk." ));
388 0 : }
389 :
390 : // Get/check shapes ( Assume that the PSF always exists in the imstore... A valid assumption in the context of mtmfs_via_cube )
391 0 : IPosition cube_shp( cube_imstore->psf()->shape() );
392 0 : IPosition mt_shp( mt_imstore->psf(0)->shape() );
393 0 : if( cube_shp[0] != mt_shp[0] || cube_shp[1] != mt_shp[1] ||cube_shp[2] != mt_shp[2] ){
394 0 : throw( AipsError("The Cube and Multi-Term images should have the same nx, ny and npol"));
395 : }
396 :
397 : // Read reference frequency
398 0 : Quantity reffreq_qa;
399 0 : Quantity::read( reffreq_qa, reffreq );
400 0 : Double refval = reffreq_qa.getValue("Hz");
401 :
402 : //cout << "ref freq : " << refval << endl;
403 :
404 : //Get the frequency list for the cube
405 0 : CoordinateSystem csys ( cube_imstore->getCSys() );
406 0 : Vector<Double> freqlist( cube_shp[3] );
407 :
408 0 : for(uInt i=0; i<csys.nCoordinates(); i++)
409 : {
410 0 : if( csys.type(i) == Coordinate::SPECTRAL )
411 : {
412 0 : SpectralCoordinate speccoord(csys.spectralCoordinate(i));
413 :
414 0 : for(Int ch=0;ch<cube_shp[3];ch++)
415 : {
416 : Double freq;
417 0 : Bool ret = speccoord.toWorld( freq, ch );
418 0 : if(ret==False) throw(AipsError("Cannot read channel frequency"));
419 0 : freqlist[ch] = freq;
420 : // cout << "freq " << ch << " is " << freq << endl;
421 : }
422 :
423 0 : }
424 : }
425 :
426 :
427 : // Reset the Taylor Sum values to zero.
428 0 : for(Int i=0;i<out_nterms;i++)
429 : {
430 0 : switch(imtype)
431 : {
432 0 : case 0:mt_imstore->psf(i)->set(0.0);break;
433 0 : case 1:mt_imstore->residual(i)->set(0.0);break;
434 0 : case 2:mt_imstore->pb(i)->set(0.0);break;
435 0 : case 3:mt_imstore->sumwt(i)->set(0.0);break;
436 : }
437 : }
438 :
439 : // Get the sumwt spectrum.
440 0 : Array<Float> lsumwt;
441 0 : cube_imstore->sumwt()->get(lsumwt, False);
442 :
443 : // Sum the weights ( or just use accumulate...)
444 0 : LatticeExprNode msum( sum( *cube_imstore->sumwt() ) );
445 0 : Float wtsum = msum.getFloat();
446 :
447 : //cerr << "perchansumwt : shape "<< lsumwt.shape() << " " << lsumwt << " sumwt "<< wtsum << endl;
448 :
449 : //Float wtsum = cube_shp[3]; // This is sum of weights, if all weights are 1.0
450 :
451 : //For each pol, do the Cube-To-Taylor calculation.
452 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
453 : {
454 0 : Vector< CountedPtr <ImageInterface<Float> > > mt_subims(out_nterms);
455 0 : for(Int i=0;i<out_nterms;i++)
456 : {
457 0 : switch(imtype)
458 : {
459 0 : case 0:use_mt=mt_imstore->psf(i);break;
460 0 : case 1:use_mt=mt_imstore->residual(i);break;
461 0 : case 2:use_mt=mt_imstore->pb(i);break;
462 0 : case 3:use_mt=mt_imstore->sumwt(i);break;
463 : }
464 0 : mt_subims[i] = mt_imstore->makeSubImage(0,1,
465 0 : 0, cube_shp[3],
466 0 : pol, cube_shp[2],
467 0 : *use_mt );
468 : }
469 :
470 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
471 : {
472 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
473 0 : chan, cube_shp[3],
474 0 : pol, cube_shp[2],
475 0 : *use_cube);
476 0 : if(imtype < 2){
477 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
478 0 : chan, cube_shp[3],
479 0 : pol, cube_shp[2],
480 0 : *(cube_imstore->pb()));
481 0 : CountedPtr<ImageInterface<Float> > tmplat = new TempImage<Float>(cube_subim->shape(), cube_subim->coordinates());
482 0 : tmplat->copyData(LatticeExpr<Float>((*pb_subim) *(*cube_subim)));
483 0 : cube_subim = tmplat;
484 :
485 0 : }
486 :
487 0 : IPosition pos(4,0,0,pol,chan);
488 :
489 0 : Double wt = (freqlist[chan] - refval) / refval;
490 :
491 0 : for(Int tt=0;tt<out_nterms;tt++)
492 : {
493 0 : Double fac = pow(wt,tt);
494 : //cerr << "BEF accum " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
495 0 : LatticeExpr<Float> eachterm = LatticeExpr<Float>( (*mt_subims[tt]) + ((fac) * (*cube_subim) * lsumwt(pos))) ;
496 0 : mt_subims[tt]->copyData(eachterm);
497 : //cerr <<" AFT accum : chan " << chan << " tt " << tt << " fac " << fac << " lsumwt " << lsumwt(pos) << " pos " << pos << " max " << max(mt_subims[tt]->get()) << endl;
498 0 : }
499 :
500 :
501 0 : }//for chan
502 :
503 :
504 : // Divide by sum of weights.
505 0 : for(Int tt=0;tt<out_nterms;tt++)
506 : {
507 : //cerr << "bef div : tt " << tt << " : " << max(mt_subims[tt]->get()) << " for imtype " << imtype << endl;
508 :
509 0 : LatticeExpr<Float> eachterm;
510 0 : if (imtype < 2) {
511 0 : eachterm = LatticeExpr<Float>( iif( (*(mt_imstore->pb(0))) > pblimit , (*mt_subims[tt]) / wtsum/(*(mt_imstore->pb(0))), 0.0));
512 : }
513 : else{
514 0 : eachterm = LatticeExpr<Float>( (*mt_subims[tt]) / wtsum ) ;
515 : }
516 0 : mt_subims[tt]->copyData(eachterm);
517 0 : mt_subims[tt]->flush();
518 : // cerr << "aft div : " << max(mt_subims[tt]->get()) << endl;
519 0 : }
520 :
521 0 : }// 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 0 : if( imtype==2 )
526 : {
527 0 : mt_imstore->removeMask( mt_imstore->pb(0) );
528 : {
529 : //MSK//
530 0 : LatticeExpr<Bool> pbmask( iif( *mt_imstore->pb(0) > fabs(pblimit) , True , False ) );
531 : //MSK//
532 0 : mt_imstore->createMask( pbmask, mt_imstore->pb(0) );
533 0 : mt_imstore->pb(0)->pixelMask().unlock();
534 0 : }
535 :
536 : }
537 :
538 0 : return True;
539 :
540 0 : }//end of func
541 :
542 :
543 0 : Bool SynthesisUtilMethods::removeFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
544 : {
545 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "removeFreqDepPB"));
546 :
547 : // Set up imstores
548 0 : CountedPtr<SIImageStore> cube_imstore;
549 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
550 :
551 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
552 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
553 :
554 : try{
555 0 : cube_imstore->residual(); // Residual Cube
556 0 : cube_imstore->pb(); // PB cube
557 0 : 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 0 : IPosition cube_shp( cube_imstore->residual()->shape() );
566 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
567 0 : 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 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
573 : {
574 :
575 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
576 0 : 0, cube_shp[3],
577 0 : pol, cube_shp[2],
578 0 : (*mt_imstore->pb(0)) );
579 :
580 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
581 0 : Float mtpbmaxval = mtpbmax.getFloat();
582 0 : if(mtpbmaxval <=0.0){os << LogIO::WARN << "pb.tt0 max is < or = zero" << LogIO::POST;}
583 :
584 :
585 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
586 : {
587 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
588 0 : chan, cube_shp[3],
589 0 : pol, cube_shp[2],
590 0 : *cube_imstore->residual() );
591 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
592 0 : chan, cube_shp[3],
593 0 : pol, cube_shp[2],
594 0 : *cube_imstore->pb() );
595 :
596 0 : LatticeExprNode pbmax( max( *pb_subim ) );
597 0 : Float pbmaxval = pbmax.getFloat();
598 0 : if( pbmaxval<=0.0 )
599 : {
600 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
601 : }
602 : else
603 : {
604 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*mt_subim)*(*(cube_subim))/(*(pb_subim)) , 0.0 ) );
605 0 : cube_subim->copyData( thepbcor );
606 0 : }// if not zero
607 :
608 0 : }//for chan
609 0 : }// for pol
610 :
611 0 : return True;
612 :
613 0 : }//end of func
614 :
615 :
616 :
617 0 : Bool SynthesisUtilMethods::applyFreqDepPB(const String& cubename, const String& mtname, const Float pblimit)
618 : {
619 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "applyFreqDepPB"));
620 :
621 : // Set up imstores
622 0 : CountedPtr<SIImageStore> cube_imstore;
623 0 : cube_imstore = CountedPtr<SIImageStore>(new SIImageStore( cubename, true, true ));
624 :
625 0 : CountedPtr<SIImageStoreMultiTerm> mt_imstore;
626 0 : mt_imstore = CountedPtr<SIImageStoreMultiTerm>(new SIImageStoreMultiTerm( mtname, 1, true, true ));
627 :
628 : try{
629 0 : cube_imstore->model(); // Model Cube
630 0 : cube_imstore->pb(); // PB cube
631 0 : 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 0 : IPosition cube_shp( cube_imstore->model()->shape() );
640 0 : IPosition mt_shp( mt_imstore->pb(0)->shape() );
641 0 : 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 0 : for(Int pol=0; pol<cube_shp[2]; pol++)
647 : {
648 0 : CountedPtr<ImageInterface<Float> > mt_subim = mt_imstore->makeSubImage(0,1,
649 0 : 0, cube_shp[3],
650 0 : pol, cube_shp[2],
651 0 : (*mt_imstore->pb(0)) );
652 :
653 0 : LatticeExprNode mtpbmax( max( *mt_subim ) );
654 0 : Float mtpbmaxval = mtpbmax.getFloat();
655 0 : 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 0 : for(Int chan=0; chan<cube_shp[3]; chan++)
661 : {
662 0 : CountedPtr<ImageInterface<Float> > cube_subim=cube_imstore->makeSubImage(0,1,
663 0 : chan, cube_shp[3],
664 0 : pol, cube_shp[2],
665 0 : *cube_imstore->model() );
666 0 : CountedPtr<ImageInterface<Float> > pb_subim=cube_imstore->makeSubImage(0,1,
667 0 : chan, cube_shp[3],
668 0 : pol, cube_shp[2],
669 0 : *cube_imstore->pb() );
670 :
671 :
672 0 : LatticeExprNode pbmax( max( *pb_subim ) );
673 0 : Float pbmaxval = pbmax.getFloat();
674 0 : if( pbmaxval<=0.0 )
675 : {
676 0 : os << LogIO::WARN << "pb max is zero for chan" << chan << LogIO::POST;
677 : }
678 : else
679 : {
680 0 : LatticeExpr<Float> thepbcor( iif( *(pb_subim) > pblimit , (*(cube_subim)) *(*(pb_subim)) / (*mt_subim) , 0.0 ) );
681 0 : cube_subim->copyData( thepbcor );
682 0 : }// if not zero
683 :
684 0 : }//for chan
685 : }// if mtpb >0
686 0 : }// for pol
687 :
688 0 : return True;
689 :
690 0 : }//end of func
691 :
692 :
693 :
694 :
695 : /***make a record of synthesisimager::weight parameters***/
696 0 : 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 0 : Record outRec;
704 0 : outRec.define("type", type);
705 0 : outRec.define("rmode", rmode);
706 0 : Record quantRec;
707 0 : QuantumHolder(noise).toRecord(quantRec);
708 0 : outRec.defineRecord("noise", quantRec);
709 0 : outRec.define("robust", robust);
710 0 : QuantumHolder(fieldofview).toRecord(quantRec);
711 0 : outRec.defineRecord("fieldofview", quantRec);
712 0 : outRec.define("npixels", npixels);
713 0 : outRec.define("multifield", multiField);
714 0 : outRec.define("usecubebriggs", useCubeBriggs);
715 0 : outRec.define("filtertype", filtertype);
716 0 : QuantumHolder(filterbmaj).toRecord(quantRec);
717 0 : outRec.defineRecord("filterbmaj", quantRec);
718 0 : QuantumHolder(filterbmin).toRecord(quantRec);
719 0 : outRec.defineRecord("filterbmin", quantRec);
720 0 : QuantumHolder(filterbpa).toRecord(quantRec);
721 0 : outRec.defineRecord("filterbpa", quantRec);
722 0 : outRec.define("fracBW", fracBW);
723 :
724 0 : return outRec;
725 0 : }
726 0 : 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 0 : QuantumHolder qh;
733 0 : String err;
734 0 : if(!inRec.isDefined("type"))
735 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
736 0 : inRec.get("type", type);
737 0 : inRec.get("rmode", rmode);
738 0 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
739 0 : throw(AipsError("Error in reading noise param"));
740 0 : noise=qh.asQuantity();
741 0 : inRec.get("robust", robust);
742 0 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
743 0 : throw(AipsError("Error in reading fieldofview param"));
744 0 : fieldofview=qh.asQuantity();
745 0 : inRec.get("npixels", npixels);
746 0 : inRec.get("multifield", multiField);
747 0 : inRec.get("usecubebriggs", useCubeBriggs);
748 0 : inRec.get("filtertype", filtertype);
749 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
750 0 : throw(AipsError("Error in reading filterbmaj param"));
751 0 : filterbmaj=qh.asQuantity();
752 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
753 0 : throw(AipsError("Error in reading filterbmin param"));
754 0 : filterbmin=qh.asQuantity();
755 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
756 0 : throw(AipsError("Error in reading filterbpa param"));
757 0 : filterbpa=qh.asQuantity();
758 0 : inRec.get("fracBW", fracBW);
759 :
760 :
761 :
762 0 : }
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 0 : 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 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "adviseChanSel"));
817 0 : if(ms==String("")){
818 0 : throw(AipsError("Need a valid MS"));
819 : }
820 0 : spw.resize();
821 0 : start.resize();
822 0 : nchan.resize();
823 : try {
824 0 : 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 0 : MeasurementSet elms(ms, TableLock(TableLock::AutoNoReadLocking), Table::Old);
856 0 : MSSelection thisSelection;
857 0 : String spsel=spwselection;
858 0 : if(spsel=="")spsel="*";
859 0 : thisSelection.setSpwExpr(spsel);
860 0 : TableExprNode exprNode=thisSelection.toTableExprNode(&elms);
861 0 : Matrix<Int> chanlist=thisSelection.getChanList();
862 0 : if(chanlist.ncolumn() <3){
863 0 : freqStart=-1.0;
864 0 : freqEnd=-1.0;
865 0 : return false;
866 : }
867 0 : Vector<Int> elspw=chanlist.column(0);
868 0 : Vector<Int> elstart=chanlist.column(1);
869 0 : Vector<Int> elnchan=Vector<Int> (chanlist.column(2)-elstart)+1;
870 0 : 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 0 : MSUtil::getFreqRangeInSpw(freqStart, freqEnd, elspw, elstart, elnchan, elms, freqframe, field_id);
892 0 : }
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 0 : return true;
912 :
913 0 : }
914 :
915 0 : 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 0 : Bool isOn = false;
923 0 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
924 0 : if (!isOn)
925 0 : 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 0 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
1397 : {
1398 0 : MVAngle mvRa=direction.getAngle().getValue()(0);
1399 0 : MVAngle mvDec=direction.getAngle().getValue()(1);
1400 0 : ostringstream oos;
1401 0 : oos << " ";
1402 0 : Int widthRA=20;
1403 0 : Int widthDec=20;
1404 0 : oos.setf(ios::left, ios::adjustfield);
1405 0 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
1406 0 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
1407 0 : oos << " "
1408 0 : << MDirection::showType(direction.getRefPtr()->getType());
1409 0 : return String(oos);
1410 0 : }
1411 :
1412 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1413 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1414 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
1415 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1416 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1417 :
1418 : // Read String from Record
1419 0 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
1420 : {
1421 0 : if( rec.isDefined( id ) )
1422 : {
1423 0 : String inval("");
1424 0 : if( rec.dataType( id )==TpString )
1425 0 : { 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 0 : if(inval.length()>0){val=inval;}
1431 0 : return String("");
1432 : }
1433 0 : else { return String(id + " must be a string\n"); }
1434 0 : }
1435 0 : else { return String("");}
1436 : }
1437 :
1438 : // Read Integer from Record
1439 0 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
1440 : {
1441 0 : if( rec.isDefined( id ) )
1442 : {
1443 0 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
1444 0 : else { return String(id + " must be an integer\n"); }
1445 : }
1446 0 : else { return String(""); }
1447 : }
1448 :
1449 : // Read Float from Record
1450 0 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
1451 : {
1452 0 : if( rec.isDefined( id ) )
1453 : {
1454 0 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
1455 0 : { rec.get( id , val ); return String(""); }
1456 0 : else { return String(id + " must be a float\n"); }
1457 : }
1458 0 : else { return String(""); }
1459 : }
1460 :
1461 : // Read Bool from Record
1462 0 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
1463 : {
1464 0 : if( rec.isDefined( id ) )
1465 : {
1466 0 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
1467 0 : else { return String(id + " must be a bool\n"); }
1468 : }
1469 0 : else{ return String(""); }
1470 : }
1471 :
1472 : // Read Vector<Int> from Record
1473 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
1474 : {
1475 0 : if( rec.isDefined( id ) )
1476 : {
1477 0 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1478 0 : { 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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1492 : {
1493 0 : if( rec.isDefined( id ) )
1494 : {
1495 0 : if( rec.dataType( id )==TpArrayFloat )
1496 : {
1497 0 : 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 0 : else if ( rec.dataType( id ) ==TpArrayDouble )
1511 : {
1512 0 : Vector<Double> vec; rec.get( id, vec );
1513 0 : val.resize(vec.nelements());
1514 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1515 0 : return String("");
1516 0 : }
1517 0 : else if ( rec.dataType( id ) ==TpArrayInt )
1518 : {
1519 0 : Vector<Int> vec; rec.get( id, vec );
1520 0 : val.resize(vec.nelements());
1521 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1522 0 : return String("");
1523 0 : }
1524 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1525 : {
1526 0 : Vector<Bool> vec; rec.get( id, vec );
1527 0 : 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 0 : }
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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1538 : {
1539 0 : if( rec.isDefined( id ) )
1540 : {
1541 0 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1542 0 : { rec.get( id , val ); return String(""); }
1543 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1544 : {
1545 0 : Vector<Bool> vec; rec.get( id, vec );
1546 0 : if( vec.nelements()==0 ){val.resize(0); return String("");}
1547 0 : else{ return String(id + " must be a vector of strings.\n"); }
1548 0 : }
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 0 : 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 0 : 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 0 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1569 : {
1570 : try
1571 : {
1572 0 : String tmpRF, tmpRA, tmpDEC;
1573 0 : std::istringstream iss(instr);
1574 0 : iss >> tmpRF >> tmpRA >> tmpDEC;
1575 0 : 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 0 : casacore::Quantity tmpQRA;
1582 0 : casacore::Quantity tmpQDEC;
1583 0 : casacore::Quantity::read(tmpQRA, tmpRA);
1584 0 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1585 :
1586 0 : 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 0 : Bool status = MDirection::getType(theRF, tmpRF);
1597 0 : if (!status) {
1598 0 : throw AipsError();
1599 : }
1600 0 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1601 0 : return String("");
1602 0 : }
1603 0 : catch(AipsError &x)
1604 : {
1605 0 : return String("Error in converting '" + instr + "' to MDirection. Need format : 'J2000 19:59:28.500 +40.44.01.50'\n");
1606 0 : }
1607 : }
1608 :
1609 : // Read Quantity from a Record string
1610 0 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1611 : {
1612 0 : if( rec.isDefined( id ) )
1613 : {
1614 0 : if( rec.dataType( id )==TpString )
1615 0 : { 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 0 : else{ return String(""); }
1619 : }
1620 :
1621 : // Read MDirection from a Record string
1622 0 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1623 : {
1624 0 : if( rec.isDefined( id ) )
1625 : {
1626 0 : if( rec.dataType( id )==TpString )
1627 0 : { 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 0 : else{ return String(""); }
1631 : }
1632 :
1633 : // Convert MDirection to String
1634 0 : String SynthesisParams::MDirectionToString(MDirection val) const
1635 : {
1636 0 : MVDirection mvpc( val.getAngle() );
1637 0 : MVAngle ra = mvpc.get()(0);
1638 0 : MVAngle dec = mvpc.get()(1);
1639 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1640 0 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1641 0 : }
1642 :
1643 : // Convert Quantity to String
1644 0 : String SynthesisParams::QuantityToString(Quantity val) const
1645 : {
1646 0 : 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 0 : ss.precision(std::numeric_limits<double>::digits10+2);
1651 0 : ss << val;
1652 0 : 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 0 : }
1658 :
1659 : // Convert Record contains Quantity or Measure quantities to String
1660 0 : String SynthesisParams::recordQMToString(const Record &rec) const
1661 : {
1662 : Double val;
1663 0 : String unit;
1664 0 : if ( rec.isDefined("m0") )
1665 : {
1666 0 : Record subrec = rec.subRecord("m0");
1667 0 : subrec.get("value",val);
1668 0 : subrec.get("unit",unit);
1669 0 : }
1670 0 : else if (rec.isDefined("value") )
1671 : {
1672 0 : rec.get("value",val);
1673 0 : rec.get("unit",unit);
1674 : }
1675 0 : return String::toString(val) + unit;
1676 0 : }
1677 :
1678 :
1679 : /////////////////////// Selection Parameters
1680 :
1681 0 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1682 : {
1683 0 : setDefaults();
1684 0 : }
1685 :
1686 0 : SynthesisParamsSelect::~SynthesisParamsSelect()
1687 : {
1688 0 : }
1689 :
1690 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1691 0 : operator=(other);
1692 0 : }
1693 :
1694 0 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1695 0 : if(this!=&other) {
1696 0 : msname=other.msname;
1697 0 : spw=other.spw;
1698 0 : freqbeg=other.freqbeg;
1699 0 : freqend=other.freqend;
1700 0 : freqframe=other.freqframe;
1701 0 : field=other.field;
1702 0 : antenna=other.antenna;
1703 0 : timestr=other.timestr;
1704 0 : scan=other.scan;
1705 0 : obs=other.obs;
1706 0 : state=other.state;
1707 0 : uvdist=other.uvdist;
1708 0 : taql=other.taql;
1709 0 : usescratch=other.usescratch;
1710 0 : readonly=other.readonly;
1711 0 : incrmodel=other.incrmodel;
1712 0 : datacolumn=other.datacolumn;
1713 :
1714 : }
1715 0 : return *this;
1716 : }
1717 :
1718 0 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1719 : {
1720 0 : setDefaults();
1721 :
1722 0 : String err("");
1723 :
1724 : try
1725 : {
1726 :
1727 0 : err += readVal( inrec, String("msname"), msname );
1728 :
1729 0 : err += readVal( inrec, String("readonly"), readonly );
1730 0 : err += readVal( inrec, String("usescratch"), usescratch );
1731 :
1732 : // override with entries from savemodel.
1733 0 : String savemodel("");
1734 0 : err += readVal( inrec, String("savemodel"), savemodel );
1735 0 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1736 0 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1737 0 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1738 :
1739 0 : err += readVal( inrec, String("incrmodel"), incrmodel );
1740 :
1741 0 : err += readVal( inrec, String("spw"), spw );
1742 :
1743 : /// -------------------------------------------------------------------------------------
1744 : // Why are these params here ? Repeats of defineimage.
1745 0 : err += readVal( inrec, String("freqbeg"), freqbeg);
1746 0 : err += readVal( inrec, String("freqend"), freqend);
1747 :
1748 0 : String freqframestr( MFrequency::showType(freqframe) );
1749 0 : err += readVal( inrec, String("outframe"), freqframestr);
1750 0 : if( ! MFrequency::getType(freqframe, freqframestr) )
1751 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1752 : /// -------------------------------------------------------------------------------------
1753 :
1754 0 : err += readVal( inrec, String("field"),field);
1755 0 : err += readVal( inrec, String("antenna"),antenna);
1756 0 : err += readVal( inrec, String("timestr"),timestr);
1757 0 : err += readVal( inrec, String("scan"),scan);
1758 0 : err += readVal( inrec, String("obs"),obs);
1759 0 : err += readVal( inrec, String("state"),state);
1760 0 : err += readVal( inrec, String("uvdist"),uvdist);
1761 0 : err += readVal( inrec, String("taql"),taql);
1762 :
1763 0 : err += readVal( inrec, String("datacolumn"),datacolumn);
1764 :
1765 0 : err += verify();
1766 :
1767 0 : }
1768 0 : catch(AipsError &x)
1769 : {
1770 0 : err = err + x.getMesg() + "\n";
1771 0 : }
1772 :
1773 0 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1774 :
1775 0 : }
1776 :
1777 0 : String SynthesisParamsSelect::verify() const
1778 : {
1779 0 : String err;
1780 : // Does the MS exist on disk.
1781 0 : Directory thems( msname );
1782 0 : if( thems.exists() )
1783 : {
1784 : // Is it readable ?
1785 0 : if( ! thems.isReadable() )
1786 0 : { err += "MS " + msname + " is not readable.\n"; }
1787 : // Depending on 'readonly', is the MS writable ?
1788 0 : 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 0 : return err;
1795 0 : }
1796 :
1797 0 : void SynthesisParamsSelect::setDefaults()
1798 : {
1799 0 : msname="";
1800 0 : spw="";
1801 0 : freqbeg="";
1802 0 : freqend="";
1803 0 : MFrequency::getType(freqframe,"LSRK");
1804 0 : field="";
1805 0 : antenna="";
1806 0 : timestr="";
1807 0 : scan="";
1808 0 : obs="";
1809 0 : state="";
1810 0 : uvdist="";
1811 0 : taql="";
1812 0 : usescratch=false;
1813 0 : readonly=true;
1814 0 : incrmodel=false;
1815 0 : datacolumn="corrected";
1816 0 : }
1817 :
1818 0 : Record SynthesisParamsSelect::toRecord()const
1819 : {
1820 0 : Record selpar;
1821 0 : selpar.define("msname",msname);
1822 0 : selpar.define("spw",spw);
1823 0 : selpar.define("freqbeg",freqbeg);
1824 0 : selpar.define("freqend",freqend);
1825 0 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1826 : //looks like fromRecord looks for outframe !
1827 0 : selpar.define("outframe", MFrequency::showType(freqframe));
1828 0 : selpar.define("field",field);
1829 0 : selpar.define("antenna",antenna);
1830 0 : selpar.define("timestr",timestr);
1831 0 : selpar.define("scan",scan);
1832 0 : selpar.define("obs",obs);
1833 0 : selpar.define("state",state);
1834 0 : selpar.define("uvdist",uvdist);
1835 0 : selpar.define("taql",taql);
1836 0 : selpar.define("usescratch",usescratch);
1837 0 : selpar.define("readonly",readonly);
1838 0 : selpar.define("incrmodel",incrmodel);
1839 0 : selpar.define("datacolumn",datacolumn);
1840 :
1841 0 : return selpar;
1842 0 : }
1843 :
1844 :
1845 : /////////////////////// Image Parameters
1846 :
1847 0 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1848 : {
1849 0 : setDefaults();
1850 0 : }
1851 :
1852 0 : SynthesisParamsImage::~SynthesisParamsImage()
1853 : {
1854 0 : }
1855 :
1856 0 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1857 0 : if(this != &other){
1858 0 : imageName=other.imageName;
1859 0 : stokes=other.stokes;
1860 0 : startModel.resize(); startModel=other.startModel;
1861 0 : imsize.resize(); imsize=other.imsize;
1862 0 : cellsize.resize(); cellsize=other.cellsize;
1863 0 : projection=other.projection;
1864 0 : useNCP=other.useNCP;
1865 0 : phaseCenter=other.phaseCenter;
1866 0 : phaseCenterFieldId=other.phaseCenterFieldId;
1867 0 : obslocation=other.obslocation;
1868 0 : pseudoi=other.pseudoi;
1869 0 : nchan=other.nchan;
1870 0 : nTaylorTerms=other.nTaylorTerms;
1871 0 : chanStart=other.chanStart;
1872 0 : chanStep=other.chanStep;
1873 0 : freqStart=other.freqStart;
1874 0 : freqStep=other.freqStep;
1875 0 : refFreq=other.refFreq;
1876 0 : velStart=other.velStart;
1877 0 : velStep=other.velStep;
1878 0 : freqFrame=other.freqFrame;
1879 0 : mFreqStart=other.mFreqStart;
1880 0 : mFreqStep=other.mFreqStep;
1881 0 : mVelStart=other.mVelStart;
1882 0 : mVelStep=other.mVelStep;
1883 0 : restFreq.resize(); restFreq=other.restFreq;
1884 0 : start=other.start;
1885 0 : step=other.step;
1886 0 : frame=other.frame;
1887 0 : veltype=other.veltype;
1888 0 : mode=other.mode;
1889 0 : reffreq=other.reffreq;
1890 0 : sysvel=other.sysvel;
1891 0 : sysvelframe=other.sysvelframe;
1892 0 : sysvelvalue=other.sysvelvalue;
1893 0 : qmframe=other.qmframe;
1894 0 : mveltype=other.mveltype;
1895 0 : tststr=other.tststr;
1896 0 : startRecord=other.startRecord;
1897 0 : stepRecord=other.stepRecord;
1898 0 : reffreqRecord=other.reffreqRecord;
1899 0 : sysvelRecord=other.sysvelRecord;
1900 0 : restfreqRecord=other.restfreqRecord;
1901 0 : csysRecord=other.csysRecord;
1902 0 : csys=other.csys;
1903 0 : imshape.resize(); imshape=other.imshape;
1904 0 : freqFrameValid=other.freqFrameValid;
1905 0 : overwrite=other.overwrite;
1906 0 : deconvolver=other.deconvolver;
1907 0 : distance=other.distance;
1908 0 : trackDir=other.trackDir;
1909 0 : trackSource=other.trackSource;
1910 0 : movingSource=other.movingSource;
1911 :
1912 :
1913 :
1914 : }
1915 :
1916 0 : return *this;
1917 :
1918 : }
1919 :
1920 0 : void SynthesisParamsImage::fromRecord(const Record &inrec)
1921 : {
1922 0 : setDefaults();
1923 0 : String err("");
1924 :
1925 : try
1926 : {
1927 :
1928 0 : err += readVal( inrec, String("imagename"), imageName);
1929 :
1930 : //// imsize
1931 0 : if( inrec.isDefined("imsize") )
1932 : {
1933 0 : DataType tp = inrec.dataType("imsize");
1934 :
1935 0 : if( tp == TpInt ) // A single integer for both dimensions.
1936 : {
1937 0 : Int npix; inrec.get("imsize", npix);
1938 0 : imsize.resize(2);
1939 0 : imsize.set( npix );
1940 : }
1941 0 : else if( tp == TpArrayInt ) // An integer array : [ nx ] or [ nx, ny ]
1942 : {
1943 0 : Vector<Int> ims;
1944 0 : inrec.get("imsize", ims);
1945 0 : if( ims.nelements()==1 ) // [ nx ]
1946 0 : {imsize.set(ims[0]); }
1947 0 : else if( ims.nelements()==2 ) // [ nx, ny ]
1948 0 : { imsize[0]=ims[0]; 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 0 : }
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 0 : if( inrec.isDefined("cell") )
1959 : {
1960 : try
1961 : {
1962 0 : DataType tp = inrec.dataType("cell");
1963 0 : if( tp == TpInt ||
1964 0 : tp == TpFloat ||
1965 : tp == TpDouble )
1966 : {
1967 0 : Double cell = inrec.asDouble("cell");
1968 0 : cellsize.set( Quantity( cell , "arcsec" ) );
1969 0 : }
1970 0 : else if ( tp == TpArrayInt ||
1971 0 : tp == TpArrayFloat ||
1972 : tp == TpArrayDouble )
1973 : {
1974 0 : Vector<Double> cells;
1975 0 : inrec.get("cell", cells);
1976 0 : if(cells.nelements()==1) // [ cellx ]
1977 0 : {cellsize.set( Quantity( cells[0], "arcsec" ) ); }
1978 0 : else if( cells.nelements()==2 ) // [ cellx, celly ]
1979 0 : { 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 0 : }
1983 0 : else if( tp == TpString )
1984 : {
1985 0 : String cell;
1986 0 : inrec.get("cell",cell);
1987 0 : Quantity qcell;
1988 0 : err += stringToQuantity( cell, qcell );
1989 0 : cellsize.set( qcell );
1990 0 : }
1991 0 : else if( tp == TpArrayString )
1992 : {
1993 0 : Array<String> cells;
1994 0 : inrec.get("cell", cells);
1995 0 : Vector<String> vcells(cells);
1996 0 : if(cells.nelements()==1) // [ cellx ]
1997 : {
1998 0 : Quantity qcell;
1999 0 : err+= stringToQuantity( vcells[0], qcell ); cellsize.set( qcell );
2000 0 : }
2001 0 : else if( cells.nelements()==2 ) // [ cellx, celly ]
2002 : {
2003 0 : err+= stringToQuantity( vcells[0], cellsize[0] );
2004 0 : 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 0 : }
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 0 : err += readVal( inrec, String("stokes"), stokes);
2021 0 : if(stokes.matches("pseudoI"))
2022 : {
2023 0 : stokes="I";
2024 0 : pseudoi=true;
2025 : }
2026 0 : else {pseudoi=false;}
2027 :
2028 : /// PseudoI
2029 :
2030 : ////nchan
2031 0 : err += readVal( inrec, String("nchan"), nchan);
2032 :
2033 : /// phaseCenter (as a string) . // Add INT support later.
2034 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
2035 0 : if( inrec.isDefined("phasecenter") )
2036 : {
2037 0 : String pcent("");
2038 0 : if( inrec.dataType("phasecenter") == TpString )
2039 : {
2040 0 : inrec.get("phasecenter",pcent);
2041 0 : if( pcent.length() > 0 ) // if it's zero length, it means 'figure out from first field in MS'.
2042 : {
2043 0 : err += readVal( inrec, String("phasecenter"), phaseCenter );
2044 0 : 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 0 : err += readVal( inrec, String("phasecenterfieldid"), phaseCenterFieldId);
2048 : }
2049 : //else { phaseCenterFieldId=0; } // Take the first field of the MS.
2050 0 : else { phaseCenterFieldId=-2; } // deal with this later in buildCoordinateSystem to assign the first selected field
2051 : }
2052 0 : if (inrec.dataType("phasecenter")==TpInt
2053 0 : || inrec.dataType("phasecenter")==TpFloat
2054 0 : || inrec.dataType("phasecenter")==TpDouble )
2055 : {
2056 : // This will override the previous setting to 0 if the phaseCenter string is zero length.
2057 0 : err += readVal( inrec, String("phasecenter"), phaseCenterFieldId );
2058 : }
2059 :
2060 0 : if( ( inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter")!=TpInt
2061 0 : && 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 0 : }
2067 :
2068 :
2069 : //// Projection
2070 0 : if( inrec.isDefined("projection") )
2071 : {
2072 0 : if( inrec.dataType("projection") == TpString )
2073 : {
2074 0 : String pstr;
2075 0 : inrec.get("projection",pstr);
2076 :
2077 : try
2078 : {
2079 0 : if( pstr.matches("NCP") )
2080 : {
2081 0 : pstr ="SIN";
2082 0 : useNCP=true;
2083 : }
2084 0 : projection=Projection::type( pstr );
2085 : }
2086 0 : catch(AipsError &x)
2087 : {
2088 0 : err += String("Invalid projection code : " + pstr );
2089 0 : }
2090 0 : }
2091 0 : else { err += "projection must be a string\n"; }
2092 : }//projection
2093 :
2094 : // Frequency frame stuff.
2095 0 : err += readVal( inrec, String("specmode"), mode);
2096 : // Alias for 'mfs' is 'cont'
2097 0 : if(mode=="cont") mode="mfs";
2098 :
2099 0 : err += readVal( inrec, String("outframe"), frame);
2100 0 : qmframe="";
2101 : // mveltype is only set when start/step is given in mdoppler
2102 0 : mveltype="";
2103 : //start
2104 0 : String startType("");
2105 0 : String widthType("");
2106 0 : if( inrec.isDefined("start") )
2107 : {
2108 0 : if( inrec.dataType("start") == TpInt )
2109 : {
2110 0 : err += readVal( inrec, String("start"), chanStart);
2111 0 : start = String::toString(chanStart);
2112 0 : startType="chan";
2113 : }
2114 0 : else if( inrec.dataType("start") == TpString )
2115 : {
2116 0 : err += readVal( inrec, String("start"), start);
2117 0 : if( start.contains("Hz") )
2118 : {
2119 0 : stringToQuantity(start,freqStart);
2120 0 : startType="freq";
2121 : }
2122 0 : else if( start.contains("m/s") )
2123 : {
2124 0 : stringToQuantity(start,velStart);
2125 0 : startType="vel";
2126 : }
2127 : }
2128 0 : 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 0 : startRecord = inrec.subRecord("start");
2134 0 : if(startRecord.isDefined("m0") )
2135 : {
2136 : //must be a measure
2137 0 : String mtype;
2138 0 : startRecord.get("type", mtype);
2139 0 : if( mtype=="frequency")
2140 : {
2141 : //mfrequency
2142 0 : startRecord.get(String("refer"), qmframe);
2143 0 : 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 0 : start = recordQMToString(startRecord);
2149 0 : stringToQuantity(start,freqStart);
2150 0 : startType="freq";
2151 : }
2152 0 : else if( mtype=="radialvelocity")
2153 : {
2154 : //mradialvelocity
2155 0 : startRecord.get(String("refer"), qmframe);
2156 0 : if ( frame!="" && frame!=qmframe)
2157 : {
2158 : // should emit warning to the logger
2159 0 : cerr<<"The frame in start:"<<qmframe<<" Override frame="<<frame<<endl;
2160 : }
2161 0 : start = recordQMToString(startRecord);
2162 0 : stringToQuantity(start,velStart);
2163 0 : startType="vel";
2164 : }
2165 0 : else if( mtype=="doppler")
2166 : {
2167 : //use veltype in mdoppler
2168 : //start = MDopToVelString(startRecord);
2169 0 : start = recordQMToString(startRecord);
2170 0 : stringToQuantity(start,velStart);
2171 0 : startRecord.get(String("refer"), mveltype);
2172 0 : mveltype.downcase();
2173 0 : startType="vel";
2174 : }
2175 0 : }
2176 : else
2177 : {
2178 0 : start = recordQMToString(startRecord);
2179 0 : if ( start.contains("Hz") )
2180 : {
2181 0 : stringToQuantity(start,freqStart);
2182 0 : startType="freq";
2183 : }
2184 0 : else if ( start.contains("m/s") )
2185 : {
2186 0 : stringToQuantity(start,velStart);
2187 0 : 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 0 : if( inrec.isDefined("width") )
2197 : {
2198 0 : if( inrec.dataType("width") == TpInt )
2199 : {
2200 0 : err += readVal( inrec, String("width"), chanStep);
2201 0 : step = String::toString(chanStep);
2202 0 : widthType="chan";
2203 : }
2204 0 : else if( inrec.dataType("width") == TpString )
2205 : {
2206 0 : err += readVal( inrec, String("width"), step);
2207 0 : if( step.contains("Hz") )
2208 : {
2209 0 : stringToQuantity(step,freqStep);
2210 0 : widthType="freq";
2211 : }
2212 0 : else if( step.contains("m/s") )
2213 : {
2214 0 : stringToQuantity(step,velStep);
2215 0 : widthType="vel";
2216 : }
2217 : }
2218 0 : 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 0 : stepRecord = inrec.subRecord("width");
2224 0 : if(stepRecord.isDefined("m0") )
2225 : {
2226 : //must be a measure
2227 0 : String mtype;
2228 0 : stepRecord.get("type", mtype);
2229 0 : 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 0 : else if( mtype=="radialvelocity")
2243 : {
2244 : //mradialvelocity
2245 0 : stepRecord.get(String("refer"), qmframe);
2246 0 : 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 0 : step = recordQMToString(stepRecord);
2252 0 : stringToQuantity(step,velStep);
2253 0 : 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 0 : }
2265 : else
2266 : {
2267 0 : step = recordQMToString(stepRecord);
2268 0 : if ( step.contains("Hz") )
2269 : {
2270 0 : stringToQuantity(step,freqStep);
2271 0 : widthType="freq";
2272 : }
2273 0 : else if ( step.contains("m/s") )
2274 : {
2275 0 : stringToQuantity(step,velStep);
2276 0 : 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 0 : 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 0 : if( inrec.isDefined("reffreq") )
2290 : {
2291 0 : if( inrec.dataType("reffreq")==TpString )
2292 : {
2293 0 : 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 0 : err += readVal( inrec, String("veltype"), veltype);
2321 0 : veltype = mveltype!=""? mveltype:veltype;
2322 : // sysvel (String, Quantity)
2323 0 : if( inrec.isDefined("sysvel") )
2324 : {
2325 0 : if( inrec.dataType("sysvel")==TpString )
2326 : {
2327 0 : 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 0 : err += readVal( inrec, String("sysvelframe"), sysvelframe);
2340 :
2341 : // rest frequencies (record or vector of Strings)
2342 0 : if( inrec.isDefined("restfreq") )
2343 : {
2344 0 : Vector<String> rfreqs(0);
2345 0 : Record restfreqSubRecord;
2346 0 : 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 0 : else if( inrec.dataType("restfreq")==TpArrayString )
2361 : {
2362 : //Vector<String> rfreqs(0);
2363 0 : err += readVal( inrec, String("restfreq"), rfreqs );
2364 : // case no restfreq is given: set to
2365 : }
2366 0 : else if( inrec.dataType("restfreq")==TpString )
2367 : {
2368 0 : rfreqs.resize(1);
2369 0 : err += readVal( inrec, String("restfreq"), rfreqs[0] );
2370 : // case no restfreq is given: set to
2371 : }
2372 0 : restFreq.resize( rfreqs.nelements() );
2373 0 : for( uInt fr=0; fr<rfreqs.nelements(); fr++)
2374 : {
2375 0 : err += stringToQuantity( rfreqs[fr], restFreq[fr] );
2376 : }
2377 0 : } // if def restfreq
2378 :
2379 : // optional - coordsys, imshape
2380 : // if exist use them. May need a consistency check with the rest of impars?
2381 0 : if( inrec.isDefined("csys") )
2382 : {
2383 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
2384 0 : if( inrec.dataType("csys")==TpRecord )
2385 : {
2386 : //csysRecord = inrec.subRecord("csys");
2387 0 : csysRecord.defineRecord("coordsys",inrec.subRecord("csys"));
2388 : }
2389 0 : if( inrec.isDefined("imshape") )
2390 : {
2391 0 : if ( inrec.dataType("imshape") == TpArrayInt )
2392 : {
2393 0 : 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 0 : String freqframestr = (frame!="" && qmframe!="")? qmframe:frame;
2404 0 : if( frame!="" && ! MFrequency::getType(freqFrame, freqframestr) )
2405 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
2406 0 : err += readVal( inrec, String("restart"), overwrite );
2407 :
2408 0 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
2409 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
2410 0 : if( inrec.isDefined("startmodel") )
2411 : {
2412 0 : if( inrec.dataType("startmodel")==TpString )
2413 : {
2414 0 : String onemodel;
2415 0 : err += readVal( inrec, String("startmodel"), onemodel );
2416 0 : if( onemodel.length()>0 )
2417 : {
2418 0 : startModel.resize(1);
2419 0 : startModel[0] = onemodel;
2420 : }
2421 0 : else {startModel.resize();}
2422 0 : }
2423 0 : else if( inrec.dataType("startmodel")==TpArrayString ||
2424 0 : inrec.dataType("startmodel")==TpArrayBool)
2425 : {
2426 0 : 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 0 : err += readVal( inrec, String("nterms"), nTaylorTerms );
2434 0 : err += readVal( inrec, String("deconvolver"), deconvolver );
2435 :
2436 : // Force nchan=1 for anything other than cube modes...
2437 0 : if(mode=="mfs") nchan=1;
2438 : //read obslocation
2439 0 : if(inrec.isDefined("obslocation_rec")){
2440 0 : String errorobs;
2441 0 : const Record obsrec=inrec.asRecord("obslocation_rec");
2442 0 : MeasureHolder mh;
2443 0 : if(!mh.fromRecord(errorobs, obsrec)){
2444 0 : err+=errorobs;
2445 : }
2446 0 : obslocation=mh.asMPosition();
2447 :
2448 0 : }
2449 :
2450 :
2451 :
2452 0 : err += verify();
2453 :
2454 0 : }
2455 0 : catch(AipsError &x)
2456 : {
2457 0 : err = err + x.getMesg() + "\n";
2458 0 : }
2459 :
2460 0 : if( err.length()>0 ) throw(AipsError("Invalid Image Parameter set : " + err));
2461 :
2462 0 : }
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 0 : String SynthesisParamsImage::verify() const
2496 : {
2497 0 : String err;
2498 :
2499 0 : if( imageName=="" ) {err += "Please supply an image name\n";}
2500 :
2501 0 : if( imsize.nelements() != 2 ){ err += "imsize must be a vector of 2 Ints\n"; }
2502 0 : if( cellsize.nelements() != 2 ) { err += "cellsize must be a vector of 2 Quantities\n"; }
2503 0 : if( cellsize[0].getValue() == 0.0 || cellsize[1].getValue() == 0.0 ) {
2504 0 : 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 0 : if( (mode=="mfs") && nchan>1 )
2512 0 : { err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";}
2513 :
2514 0 : if( ! stokes.matches("I") && ! stokes.matches("Q") &&
2515 0 : ! stokes.matches("U") && ! stokes.matches("V") &&
2516 0 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2517 0 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2518 0 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2519 0 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2520 0 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2521 0 : ! stokes.matches("IQU") && ! stokes.matches("IUV") &&
2522 0 : ! 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 0 : if( startModel.nelements()>0 )
2529 : {
2530 0 : if( deconvolver!="mtmfs" ) {
2531 :
2532 0 : if( startModel.nelements()!=1 ){err += String("Only one startmodel image is allowed.\n");}
2533 : else
2534 : {
2535 0 : File fp( imageName+String(".model") );
2536 0 : if( fp.exists() ) err += "Model " + imageName+".model exists, but a starting model of " + startModel[0] + " is also being requested. Please either reset startmodel='' to use what already exists, or delete " + imageName + ".model so that it uses the new model specified in startmodel.";
2537 0 : }
2538 : }
2539 : else {// mtmfs
2540 0 : File fp( imageName+String(".model.tt0") );
2541 0 : 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 0 : }
2547 :
2548 : // Check that startmodel exists on disk !
2549 0 : for(uInt ss=0;ss<startModel.nelements();ss++)
2550 : {
2551 0 : File fp( startModel[ss] );
2552 0 : if( ! fp.exists() ) {err += "Startmodel " + startModel[ss] + " cannot be found on disk.";}
2553 0 : }
2554 :
2555 : }
2556 :
2557 :
2558 : /// Check imsize for efficiency.
2559 0 : Int imxnew = SynthesisUtilMethods::getOptimumSize( imsize[0] );
2560 0 : Int imynew = SynthesisUtilMethods::getOptimumSize( imsize[1] );
2561 :
2562 0 : if( imxnew != imsize[0] || imynew != imsize[1] )
2563 : {
2564 0 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2565 0 : 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 0 : 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 0 : }
2569 :
2570 0 : 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 0 : void SynthesisParamsImage::setDefaults()
2584 : {
2585 : // Image definition parameters
2586 0 : imageName = String("");
2587 0 : imsize.resize(2); imsize.set(100);
2588 0 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2589 0 : stokes="I";
2590 0 : phaseCenter=MDirection();
2591 0 : phaseCenterFieldId=-1;
2592 0 : projection=Projection::SIN;
2593 0 : useNCP=false;
2594 0 : startModel=Vector<String>(0);
2595 0 : freqFrameValid=True;
2596 0 : overwrite=false;
2597 : // PseudoI
2598 0 : pseudoi=false;
2599 :
2600 : // Spectral coordinates
2601 0 : nchan=1;
2602 0 : mode="mfs";
2603 0 : start="";
2604 0 : step="";
2605 0 : chanStart=0;
2606 0 : 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 0 : freqStart=Quantity(0,"");
2612 0 : freqStep=Quantity(0,"");
2613 0 : velStart=Quantity(0,"");
2614 0 : velStep=Quantity(0,"");
2615 0 : veltype=String("radio");
2616 0 : restFreq.resize(0);
2617 0 : refFreq = Quantity(0,"Hz");
2618 0 : frame = "";
2619 0 : freqFrame=MFrequency::LSRK;
2620 0 : sysvel="";
2621 0 : sysvelframe="";
2622 0 : sysvelvalue=Quantity(0.0,"m/s");
2623 0 : nTaylorTerms=1;
2624 0 : deconvolver="hogbom";
2625 : ///csysRecord=Record();
2626 : //
2627 :
2628 :
2629 0 : }
2630 :
2631 0 : Record SynthesisParamsImage::toRecord() const
2632 : {
2633 0 : Record impar;
2634 0 : impar.define("imagename", imageName);
2635 0 : impar.define("imsize", imsize);
2636 0 : Vector<String> cells(2);
2637 0 : cells[0] = QuantityToString( cellsize[0] );
2638 0 : cells[1] = QuantityToString( cellsize[1] );
2639 0 : impar.define("cell", cells );
2640 0 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2641 0 : else{impar.define("stokes", stokes);}
2642 0 : impar.define("nchan", nchan);
2643 0 : impar.define("nterms", nTaylorTerms);
2644 0 : impar.define("deconvolver",deconvolver);
2645 0 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2646 0 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2647 0 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2648 :
2649 0 : impar.define("specmode", mode );
2650 : // start and step can be one of these types
2651 0 : if( start!="" )
2652 : {
2653 0 : if( !start.contains("Hz") && !start.contains("m/s") &&
2654 0 : String::toInt(start) == chanStart )
2655 : {
2656 0 : impar.define("start",chanStart);
2657 : }
2658 0 : else if( startRecord.nfields() > 0 )
2659 : {
2660 0 : impar.defineRecord("start", startRecord );
2661 : }
2662 : else
2663 : {
2664 0 : impar.define("start",start);
2665 : }
2666 : }
2667 : else {
2668 0 : impar.define("start", start );
2669 : }
2670 0 : if( step!="" )
2671 : {
2672 0 : if( !step.contains("Hz") && !step.contains("m/s") &&
2673 0 : String::toInt(step) == chanStep )
2674 : {
2675 0 : impar.define("width", chanStep);
2676 : }
2677 0 : else if( stepRecord.nfields() > 0 )
2678 : {
2679 0 : impar.defineRecord("width",stepRecord);
2680 : }
2681 : else
2682 : {
2683 0 : impar.define("width",step);
2684 : }
2685 : }
2686 : else
2687 : {
2688 0 : 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 0 : impar.define("veltype", veltype);
2697 0 : if (restfreqRecord.nfields() != 0 )
2698 : {
2699 0 : impar.defineRecord("restfreq", restfreqRecord);
2700 : }
2701 : else
2702 : {
2703 0 : Vector<String> rfs( restFreq.nelements() );
2704 0 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2705 0 : impar.define("restfreq", rfs);
2706 0 : }
2707 : //impar.define("reffreq", QuantityToString(refFreq));
2708 : //reffreq
2709 0 : if( reffreqRecord.nfields() != 0 )
2710 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2711 : else
2712 0 : { impar.define("reffreq",reffreq); }
2713 : //impar.define("reffreq", reffreq );
2714 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2715 0 : impar.define("outframe", frame );
2716 : //sysvel
2717 0 : if( sysvelRecord.nfields() != 0 )
2718 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2719 : else
2720 0 : { impar.define("sysvel", sysvel );}
2721 0 : impar.define("sysvelframe", sysvelframe );
2722 :
2723 0 : impar.define("restart",overwrite );
2724 0 : impar.define("freqframevalid", freqFrameValid);
2725 0 : impar.define("startmodel", startModel );
2726 :
2727 0 : if( csysRecord.isDefined("coordsys") )
2728 : {
2729 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2730 0 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2731 0 : impar.define("imshape", imshape);
2732 : }
2733 : // else cout << " NO CSYS INFO to write to output record " << endl;
2734 : ///Now save obslocation
2735 0 : Record tmprec;
2736 0 : String err;
2737 0 : MeasureHolder mh(obslocation);
2738 0 : if(mh.toRecord(err, tmprec)){
2739 0 : impar.defineRecord("obslocation_rec", tmprec);
2740 : }
2741 : else{
2742 0 : throw(AipsError("failed to save obslocation to record"));
2743 : }
2744 0 : return impar;
2745 0 : }
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 0 : 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 0 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2765 0 : vi2.originChunks();
2766 0 : vi2.origin();
2767 : /// This version uses the new vi2/vb2
2768 : // get the first ms for multiple MSes
2769 : //MeasurementSet msobj=vi2.ms();
2770 0 : Int fld=vb->fieldId()(0);
2771 :
2772 : //handling first ms only
2773 0 : Double gfreqmax=-1.0;
2774 0 : Double gdatafend=-1.0;
2775 0 : Double gdatafstart=1e14;
2776 0 : Double gfreqmin=1e14;
2777 0 : Vector<Int> spwids0;
2778 0 : Int j=0;
2779 0 : Int minfmsid=0;
2780 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2781 0 : Double imStartFreq=getCubeImageStartFreq();
2782 0 : std::vector<Int> sourceMsWithStartFreq;
2783 :
2784 :
2785 0 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2786 : //auto forMS0=chansel.find(0);
2787 0 : map<Int, Vector<Int> > spwsels=forMS0->second;
2788 0 : Int nspws=spwsels.size();
2789 0 : Vector<Int> spwids(nspws);
2790 0 : Vector<Int> nChannels(nspws);
2791 0 : Vector<Int> firstChannels(nspws);
2792 : //Vector<Int> channelIncrement(nspws);
2793 :
2794 0 : Int k=0;
2795 0 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2796 0 : spwids[k]=it->first;
2797 0 : nChannels[k]=(it->second)[0];
2798 0 : firstChannels[k]=(it->second)[1];
2799 : }
2800 0 : if(j==0) {
2801 0 : spwids0.resize();
2802 0 : 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 0 : Double freqmin=0, freqmax=0;
2814 0 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2815 :
2816 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2817 0 : 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 0 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2823 : //cerr << "after " << datafstart << " " << datafend << endl;
2824 0 : if((datafstart > datafend) || !status)
2825 0 : throw(AipsError("spw selection failed"));
2826 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2827 :
2828 0 : if (mode=="cubedata") {
2829 0 : freqmin = datafstart;
2830 0 : freqmax = datafend;
2831 : }
2832 0 : else if(mode == "cubesource"){
2833 0 : if(!trackSource){
2834 0 : throw(AipsError("Cannot be in cubesource without tracking a moving source"));
2835 : }
2836 0 : String ephemtab(movingSource);
2837 0 : if(movingSource=="TRACKFIELD"){
2838 0 : Int fieldID=MSColumns(*mss[j]).fieldId()(0);
2839 0 : ephemtab=Path(MSColumns(*mss[j]).field().ephemPath(fieldID)).absoluteName();
2840 : }
2841 0 : MEpoch refep=MSColumns(*mss[j]).timeMeas()(0);
2842 0 : Quantity refsysvel;
2843 0 : MSUtil::getFreqRangeAndRefFreqShift(freqmin,freqmax,refsysvel, refep, spwids,firstChannels, nChannels, *mss[j], ephemtab, trackDir, true);
2844 0 : if(j==0)
2845 0 : 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 0 : }
2853 : else {
2854 : //VisBufferUtil::getFreqRange(freqmin,freqmax, vi2, freqFrameValid? freqFrame:MFrequency::REST );
2855 : //cerr << "before " << freqmin << " " << freqmax << endl;
2856 0 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2857 0 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2858 : //cerr << "after " << freqmin << " " << freqmax << endl;
2859 : }
2860 :
2861 :
2862 :
2863 :
2864 0 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2865 0 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2866 0 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2867 0 : 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 0 : if(imStartFreq > 0.0 && imStartFreq >= freqmin && imStartFreq <= freqmax){
2874 0 : if(mode != "cubesource"){
2875 0 : minfmsid=j;
2876 0 : spwids0.resize();
2877 0 : spwids0=spwids;
2878 0 : vi2.originChunks();
2879 0 : vi2.origin();
2880 0 : while(vb->msId() != j && vi2.moreChunks() ){
2881 0 : vi2.nextChunk();
2882 0 : vi2.origin();
2883 : }
2884 0 : fld=vb->fieldId()(0);
2885 :
2886 : }
2887 : else{
2888 0 : sourceMsWithStartFreq.push_back(j);
2889 : }
2890 : }
2891 :
2892 0 : }
2893 0 : 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 0 : MeasurementSet msobj = *mss[minfmsid];
2900 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2901 0 : return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2902 0 : }
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 0 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2935 : MeasurementSet& msobj,
2936 : Vector<Int> spwids, Int fld,
2937 : Double freqmin, Double freqmax,
2938 : Double datafstart, Double datafend )
2939 : {
2940 0 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2941 :
2942 0 : CoordinateSystem csys;
2943 0 : if( csysRecord.nfields()!=0 )
2944 : {
2945 : //use cysRecord
2946 0 : Record subRec1;
2947 : // cout<<"USE THE EXISTING CSYS +++++++++++++++++"<<endl;
2948 0 : CoordinateSystem *csysptr = CoordinateSystem::restore(csysRecord,"coordsys");
2949 : //csys = *csysptr;
2950 : //CoordinateSystem csys(*csysptr);
2951 0 : csys = *csysptr;
2952 :
2953 0 : }
2954 : else {
2955 0 : MSColumns msc(msobj);
2956 0 : String telescop = msc.observation().telescopeName()(0);
2957 0 : MEpoch obsEpoch = msc.timeMeas()(0);
2958 0 : MPosition obsPosition;
2959 0 : 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 0 : MDirection phaseCenterToUse = phaseCenter;
2971 :
2972 0 : if( phaseCenterFieldId != -1 )
2973 : {
2974 0 : MSFieldColumns msfield(msobj.field());
2975 0 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2976 : {
2977 0 : if(trackSource){
2978 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2979 : }
2980 : else{
2981 0 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2982 : }
2983 : }
2984 : else
2985 : {
2986 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2987 : }
2988 0 : }
2989 : // Setup Phase center if it is specified only by field id.
2990 :
2991 : /////////////////// Direction Coordinates
2992 0 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2993 : // Normalize correctly
2994 0 : MVAngle ra=mvPhaseCenter.get()(0);
2995 0 : ra(0.0);
2996 :
2997 0 : MVAngle dec=mvPhaseCenter.get()(1);
2998 0 : Vector<Double> refCoord(2);
2999 0 : refCoord(0)=ra.get().getValue();
3000 0 : refCoord(1)=dec;
3001 0 : Vector<Double> refPixel(2);
3002 0 : refPixel(0) = Double(imsize[0]/2);
3003 0 : refPixel(1) = Double(imsize[1]/2);
3004 :
3005 0 : Vector<Double> deltas(2);
3006 0 : deltas(0)=-1* cellsize[0].get("rad").getValue();
3007 0 : deltas(1)=cellsize[1].get("rad").getValue();
3008 0 : Matrix<Double> xform(2,2);
3009 0 : xform=0.0;xform.diagonal()=1.0;
3010 :
3011 0 : Vector<Double> projparams(2);
3012 0 : projparams = 0.0;
3013 0 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
3014 0 : Projection projTo( projection.type(), projparams );
3015 :
3016 : DirectionCoordinate
3017 0 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
3018 : // projection,
3019 : projTo,
3020 0 : refCoord(0), refCoord(1),
3021 0 : deltas(0), deltas(1),
3022 : xform,
3023 0 : 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 0 : obslocation=obsPosition;
3031 0 : ObsInfo myobsinfo;
3032 0 : myobsinfo.setTelescope(telescop);
3033 0 : myobsinfo.setPointingCenter(mvPhaseCenter);
3034 0 : myobsinfo.setObsDate(obsEpoch);
3035 0 : 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 0 : if(spwids.nelements()==0)
3051 : {
3052 0 : Int nspw=msc.spectralWindow().nrow();
3053 0 : spwids.resize(nspw);
3054 0 : indgen(spwids);
3055 : }
3056 0 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
3057 0 : Vector<Double> dataChanFreq, dataChanWidth;
3058 0 : std::vector<std::vector<Int> > averageWhichChan;
3059 0 : std::vector<std::vector<Int> > averageWhichSPW;
3060 0 : std::vector<std::vector<Double> > averageChanFrac;
3061 :
3062 0 : if(spwids.nelements()==1)
3063 : {
3064 0 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
3065 0 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
3066 : }
3067 : else
3068 : {
3069 0 : 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 0 : Double minDataFreq = min(dataChanFreq);
3076 0 : if(start=="" && minDataFreq < datafstart ) {
3077 : // limit data chan freq vector for default start case with channel selection
3078 : Int chanStart, chanEnd;
3079 0 : Int lochan = 0;
3080 0 : Int nDataChan = dataChanFreq.nelements();
3081 0 : Int hichan = nDataChan-1;
3082 : Double diff_fmin, diff_fmax;
3083 0 : Bool ascending = dataChanFreq[nDataChan-1] - dataChanFreq[0] > 0;
3084 0 : for(Int ichan = 0; ichan < nDataChan; ichan++)
3085 : {
3086 0 : diff_fmin = dataChanFreq[ichan] - datafstart;
3087 0 : diff_fmax = datafend - dataChanFreq[ichan];
3088 : // freqmin and freqmax should corresponds to the channel edges
3089 0 : if(ascending)
3090 : {
3091 :
3092 0 : if( diff_fmin > 0 && diff_fmin <= dataChanWidth[ichan]/2. )
3093 : {
3094 0 : lochan = ichan;
3095 : }
3096 0 : else if(diff_fmax > 0 && diff_fmax <= dataChanWidth[ichan]/2. )
3097 : {
3098 0 : 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 0 : chanStart = lochan;
3114 0 : chanEnd = hichan;
3115 0 : if (lochan > hichan)
3116 : {
3117 0 : chanStart=hichan;
3118 0 : chanEnd=lochan;
3119 : }
3120 0 : Vector<Double> tempChanFreq = dataChanFreq(Slice(chanStart,chanEnd-chanStart+1,1));
3121 0 : Vector<Double> tempChanWidth = dataChanWidth(Slice(chanStart,chanEnd-chanStart+1,1));
3122 0 : dataChanFreq.resize(tempChanFreq.nelements());
3123 0 : dataChanWidth.resize(tempChanWidth.nelements());
3124 0 : dataChanFreq = tempChanFreq;
3125 0 : dataChanWidth = tempChanWidth;
3126 0 : }
3127 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3128 0 : String cubemode;
3129 0 : if ( qrestfreq.getValue("Hz")==0 )
3130 : {
3131 0 : MSDopplerUtil msdoppler(msobj);
3132 0 : Vector<Double> restfreqvec;
3133 0 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
3134 0 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
3135 0 : if ( qrestfreq.getValue("Hz")==0 and mode!="mfs" )
3136 : {
3137 0 : cubemode = findSpecMode(mode);
3138 0 : 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 0 : Double provisional_restfreq = (datafend+datafstart)/2.0;
3143 0 : 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 0 : << 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 0 : }
3153 : Double refPix;
3154 0 : Vector<Double> chanFreq;
3155 0 : Vector<Double> chanFreqStep;
3156 0 : String specmode;
3157 :
3158 0 : if(mode=="cubesource"){
3159 0 : MDoppler mdop(sysvelvalue, MDoppler::RELATIVISTIC);
3160 0 : dataChanFreq=mdop.shiftFrequency(dataChanFreq);
3161 0 : dataChanWidth=mdop.shiftFrequency(dataChanWidth);
3162 0 : 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 0 : }
3167 :
3168 0 : 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 0 : Bool nonLinearFreq(false);
3174 0 : String veltype_p=veltype;
3175 0 : veltype_p.upcase();
3176 0 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
3177 0 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
3178 : {
3179 0 : nonLinearFreq= true;
3180 : }
3181 :
3182 0 : SpectralCoordinate mySpectral;
3183 : Double stepf;
3184 0 : if(!nonLinearFreq)
3185 : //TODO: velocity mode default start case (use last channels?)
3186 : {
3187 0 : Double startf=chanFreq[0];
3188 : //Double stepf=chanFreqStep[0];
3189 0 : if(chanFreq.nelements()==1)
3190 : {
3191 0 : stepf=chanFreqStep[0];
3192 : }
3193 : else
3194 : {
3195 0 : stepf=chanFreq[1]-chanFreq[0];
3196 : }
3197 0 : Double restf=qrestfreq.getValue("Hz");
3198 : //stepf=9e8;
3199 0 : 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 0 : if(mode=="cubedata")
3203 : {
3204 : // mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
3205 0 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST?
3206 : MFrequency::REST : MFrequency::Undefined,
3207 0 : startf, stepf, refPix, restf);
3208 : }
3209 0 : 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 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
3215 0 : startf, stepf, refPix, restf);
3216 : }
3217 : else
3218 : {
3219 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3220 0 : startf, stepf, refPix, restf);
3221 : }
3222 : }
3223 : else
3224 : { // nonlinear freq coords - use tabular setting
3225 : // once NOFRAME is implemented do this
3226 0 : 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 0 : else if (mode=="cubesource")
3234 : {
3235 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
3236 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
3237 : }
3238 : else
3239 : {
3240 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
3241 0 : 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 0 : uInt nrestfreq = restFreq.nelements();
3250 0 : 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 0 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
3268 0 : if(whichStokes.nelements()==0)
3269 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
3270 0 : StokesCoordinate myStokes(whichStokes);
3271 :
3272 : ////////////////// Build Full coordinate system.
3273 :
3274 : //CoordinateSystem csys;
3275 0 : csys.addCoordinate(myRaDec);
3276 0 : csys.addCoordinate(myStokes);
3277 0 : csys.addCoordinate(mySpectral);
3278 0 : csys.setObsInfo(myobsinfo);
3279 :
3280 : //store back csys to impars record
3281 : //cerr<<"save csys to csysRecord..."<<endl;
3282 0 : if(csysRecord.isDefined("coordsys"))
3283 0 : csysRecord.removeField("coordsys");
3284 0 : csys.save(csysRecord,"coordsys");
3285 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
3286 : // imshape
3287 0 : imshape.resize(4);
3288 0 : imshape[0] = imsize[0];
3289 0 : imshape[1] = imsize[0];
3290 0 : imshape[2] = whichStokes.nelements();
3291 0 : 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 0 : } // end of else when coordsys record is not defined...
3325 :
3326 : // cout << " ----- ----- ------ ------ CSYS WORLD AXIS UNITS : " << csys.worldAxisUnits() << endl;
3327 :
3328 0 : return csys;
3329 0 : }
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 0 : MDirection SynthesisParamsImage::getMovingSourceDir(const MeasurementSet& ms, const MEpoch& refEp, const MPosition& obsposition, const MDirection::Types outframe){
3581 0 : MDirection outdir;
3582 0 : String ephemtab(movingSource);
3583 0 : if(movingSource=="TRACKFIELD"){
3584 0 : Int fieldID=MSColumns(ms).fieldId()(0);
3585 0 : ephemtab=Path(MSColumns(ms).field().ephemPath(fieldID)).absoluteName();
3586 : }
3587 0 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
3588 0 : 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 0 : MeasFrame mframe(refEp, obsposition);
3591 0 : MDirection::Ref outref1(MDirection::AZEL, mframe);
3592 0 : MDirection::Ref outref(outframe, mframe);
3593 0 : 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 0 : if (Table::isReadable(ephemtab)){
3599 0 : MeasComet mcomet(Path(ephemtab).absoluteName());
3600 0 : mframe.set(mcomet);
3601 0 : tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
3602 0 : }
3603 0 : else if (planetType >=MDirection::MERCURY && planetType <MDirection::COMET){
3604 0 : tmpazel=MDirection::Convert(trackDir, outref1)();
3605 : }
3606 0 : outdir=MDirection::Convert(tmpazel, outref)();
3607 :
3608 0 : return outdir;
3609 0 : }
3610 :
3611 0 : 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 0 : String inStart, inStep;
3622 0 : specmode = findSpecMode(mode);
3623 0 : String freqframe;
3624 0 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3625 0 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3626 :
3627 0 : refPix=0.0;
3628 0 : Bool descendingfreq(false);
3629 0 : Bool descendingoutfreq(false);
3630 :
3631 0 : if( mode.contains("cube") )
3632 : {
3633 0 : String restfreq=QuantityToString(qrestfreq);
3634 : // use frame from input start or width in MFreaquency or MRadialVelocity
3635 0 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3636 : // emit warning here if qmframe is used
3637 : //
3638 0 : inStart = start;
3639 0 : inStep = step;
3640 0 : if( specmode=="channel" )
3641 : {
3642 0 : inStart = String::toString(chanStart);
3643 0 : inStep = String::toString(chanStep);
3644 : // negative step -> descending channel indices
3645 0 : if (inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3646 : // input frame is the data frame
3647 : //freqframe = MFrequency::showType(dataFrame);
3648 : }
3649 0 : 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 0 : if(inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3660 : }
3661 0 : 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 0 : if (!inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3673 : }
3674 :
3675 0 : if (inStep=='0') inStep="";
3676 :
3677 0 : MRadialVelocity mSysVel;
3678 0 : Quantity qVel;
3679 : MRadialVelocity::Types mRef;
3680 0 : if(mode!="cubesource")
3681 : {
3682 :
3683 :
3684 0 : 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 0 : freqframe=MFrequency::showType(dataFrame);
3694 0 : 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 0 : { mSysVel=MRadialVelocity();}
3701 : }
3702 : // cubedata mode: input start, step are those of the input data frame
3703 0 : if ( mode=="cubedata" )
3704 : {
3705 0 : freqframe=MFrequency::showType(dataFrame);
3706 0 : 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 0 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3717 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3718 0 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3719 0 : << LogIO::POST;
3720 0 : ostringstream ostr;
3721 0 : ostr << " phaseCenter='" << phaseCenter;
3722 0 : os << String(ostr)<<"' ";
3723 :
3724 : Double dummy; // dummy variable - weightScale is not used here
3725 0 : 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 0 : veltype,
3742 : verbose,
3743 : mSysVel
3744 : );
3745 :
3746 0 : if( nchan==-1 )
3747 : {
3748 0 : nchan = chanFreq.nelements();
3749 0 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3750 : }
3751 :
3752 0 : if (!rst) {
3753 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
3754 0 : << LogIO::EXCEPTION;
3755 0 : return false;
3756 : }
3757 : os << LogIO::DEBUG1
3758 0 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3759 0 : << LogIO::POST;
3760 :
3761 0 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3762 0 : descendingoutfreq = true;
3763 : }
3764 :
3765 : //if (descendingfreq && !descendingoutfreq) {
3766 0 : if ((specmode=="channel" && descendingfreq==1)
3767 0 : || (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 0 : std::vector<Double> stlchanfreq;
3773 0 : chanFreq.tovector(stlchanfreq);
3774 0 : std::reverse(stlchanfreq.begin(),stlchanfreq.end());
3775 0 : chanFreq=Vector<Double>(stlchanfreq);
3776 0 : chanFreqStep=-chanFreqStep;
3777 0 : }
3778 0 : }
3779 0 : else if ( mode=="mfs" ) {
3780 0 : chanFreq.resize(1);
3781 0 : chanFreqStep.resize(1);
3782 : //chanFreqStep[0] = freqmax - freqmin;
3783 0 : Double freqmean = (freqmin + freqmax)/2;
3784 0 : if (refFreq.getValue("Hz")==0) {
3785 0 : chanFreq[0] = freqmean;
3786 0 : refPix = 0.0;
3787 0 : chanFreqStep[0] = freqmax - freqmin;
3788 : }
3789 : else {
3790 0 : chanFreq[0] = refFreq.getValue("Hz");
3791 : // Set the new reffreq to be the refPix (CAS-9518)
3792 0 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3793 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3794 0 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3795 : }
3796 :
3797 0 : if( nchan==-1 ) nchan=1;
3798 0 : if( qrestfreq.getValue("Hz")==0.0 ) {
3799 0 : restFreq.resize(1);
3800 0 : 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 0 : return true;
3810 :
3811 0 : }//getImFreq
3812 : /////////////////////////
3813 0 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3814 0 : Double inStartFreq=-1.0;
3815 0 : String checkspecmode("");
3816 0 : if(mode.contains("cube")) {
3817 0 : checkspecmode = findSpecMode(mode);
3818 : }
3819 0 : if(checkspecmode!="") {
3820 0 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3821 0 : if(checkspecmode=="channel") {
3822 0 : inStartFreq=-1.0;
3823 : }
3824 : else {
3825 0 : if(checkspecmode=="frequency") {
3826 0 : inStartFreq = freqStart.get("Hz").getValue();
3827 : }
3828 0 : else if(checkspecmode=="velocity") {
3829 : MDoppler::Types DopType;
3830 0 : MDoppler::getType(DopType, veltype);
3831 0 : MDoppler mdop(velStart,DopType);
3832 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3833 0 : inStartFreq = MFrequency::fromDoppler(mdop, qrestfreq.getValue(Unit("Hz")), mfreqframe).getValue();
3834 0 : }
3835 : }
3836 : }
3837 :
3838 0 : return inStartFreq;
3839 :
3840 0 : }
3841 :
3842 0 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3843 : {
3844 0 : String specmode;
3845 0 : specmode="channel";
3846 0 : 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 0 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3852 0 : !( velStep.getValue()==0.0)) {
3853 0 : specmode="velocity";
3854 : }
3855 0 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3856 0 : !(freqStep.getValue()==0.0)) {
3857 0 : specmode="frequency";
3858 : }
3859 : }
3860 0 : return specmode;
3861 0 : }
3862 :
3863 :
3864 0 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3865 : {
3866 0 : Vector<Int> whichStokes(0);
3867 0 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3868 0 : stokes=="RR" ||stokes=="LL" ||
3869 0 : stokes=="XX" || stokes=="YY" ) {
3870 0 : whichStokes.resize(1);
3871 0 : whichStokes(0)=Stokes::type(stokes);
3872 : }
3873 0 : else if(stokes=="IV" || stokes=="IQ" ||
3874 0 : stokes=="RRLL" || stokes=="XXYY" ||
3875 0 : stokes=="QU" || stokes=="UV"){
3876 0 : whichStokes.resize(2);
3877 :
3878 0 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3879 0 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3880 0 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3881 0 : else if(stokes=="XXYY"){whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
3882 0 : 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 0 : 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 0 : else if(stokes=="IQUV"){
3895 0 : whichStokes.resize(4);
3896 0 : whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;
3897 0 : whichStokes(2)=Stokes::U; whichStokes(3)=Stokes::V;
3898 : }
3899 :
3900 0 : return whichStokes;
3901 0 : }// decidenpolplanes
3902 :
3903 0 : IPosition SynthesisParamsImage::shp() const
3904 : {
3905 0 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3906 :
3907 0 : if( imsize[0]<=0 || imsize[1]<=0 || nStokes<=0 || nchan<=0 )
3908 : {
3909 0 : 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 0 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3913 : }
3914 :
3915 0 : Record SynthesisParamsImage::getcsys() const
3916 : {
3917 0 : 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 0 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3974 : {
3975 0 : setDefaults();
3976 0 : }
3977 :
3978 0 : SynthesisParamsGrid::~SynthesisParamsGrid()
3979 : {
3980 0 : }
3981 :
3982 :
3983 0 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3984 : {
3985 0 : setDefaults();
3986 :
3987 0 : String err("");
3988 :
3989 : try {
3990 0 : err += readVal( inrec, String("imagename"), imageName );
3991 :
3992 : // FTMachine parameters
3993 0 : err += readVal( inrec, String("gridder"), gridder );
3994 0 : err += readVal( inrec, String("padding"), padding );
3995 0 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3996 0 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3997 0 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3998 0 : err += readVal( inrec, String("convfunc"), convFunc );
3999 :
4000 0 : err += readVal( inrec, String("vptable"), vpTable );
4001 :
4002 :
4003 : // convert 'gridder' to 'ftmachine' and 'mtype'
4004 0 : ftmachine = "gridft";
4005 0 : mType = "default";
4006 0 : if (gridder=="ft" || gridder=="gridft" || gridder=="standard") {
4007 0 : ftmachine = "gridft";
4008 : }
4009 :
4010 0 : if ( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) &&
4011 0 : (wprojplanes>1 || wprojplanes==-1) ) {
4012 0 : ftmachine = "wprojectft";
4013 : }
4014 : //facetting alone use gridft
4015 0 : else if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes==1))
4016 0 : {ftmachine=="gridft";}
4017 :
4018 0 : if (gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" ) {
4019 0 : ftmachine = "mosaicft";
4020 : }
4021 :
4022 0 : if (gridder=="imagemosaic") {
4023 0 : mType = "imagemosaic";
4024 0 : if (wprojplanes>1 || wprojplanes==-1) {
4025 0 : ftmachine = "wprojectft";
4026 : }
4027 : }
4028 :
4029 0 : if (gridder=="awproject" || gridder=="awprojectft" || gridder=="awp") {
4030 0 : ftmachine = "awprojectft";
4031 : }
4032 :
4033 0 : if (gridder=="singledish") {
4034 0 : ftmachine = "sd";
4035 : }
4036 :
4037 0 : String deconvolver;
4038 0 : err += readVal( inrec, String("deconvolver"), deconvolver );
4039 0 : if (deconvolver=="mtmfs") {
4040 0 : mType = "multiterm"; // Takes precedence over imagemosaic
4041 : }
4042 :
4043 : // facets
4044 0 : err += readVal( inrec, String("facets"), facets );
4045 : // chanchunks
4046 0 : err += readVal( inrec, String("chanchunks"), chanchunks );
4047 :
4048 : // Spectral interpolation
4049 0 : err += readVal( inrec, String("interpolation"), interpolation ); // not used in SI yet...
4050 : // Track moving source ?
4051 0 : err += readVal( inrec, String("distance"), distance );
4052 0 : err += readVal( inrec, String("tracksource"), trackSource );
4053 0 : err += readVal( inrec, String("trackdir"), trackDir );
4054 :
4055 : // The extra params for WB-AWP
4056 0 : err += readVal( inrec, String("aterm"), aTermOn );
4057 0 : err += readVal( inrec, String("psterm"), psTermOn );
4058 0 : err += readVal( inrec, String("mterm"), mTermOn );
4059 0 : err += readVal( inrec, String("wbawp"), wbAWP );
4060 0 : err += readVal( inrec, String("cfcache"), cfCache );
4061 0 : err += readVal( inrec, String("usepointing"), usePointing );
4062 0 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
4063 0 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
4064 0 : err += readVal( inrec, String("conjbeams"), conjBeams );
4065 0 : err += readVal( inrec, String("computepastep"), computePAStep );
4066 0 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
4067 :
4068 : // The extra params for single-dish
4069 0 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
4070 0 : err += readVal( inrec, String("convertfirst"), convertFirst );
4071 0 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
4072 0 : err += readVal( inrec, String("convsupport"), convSupport );
4073 0 : err += readVal( inrec, String("truncate"), truncateSize );
4074 0 : err += readVal( inrec, String("gwidth"), gwidth );
4075 0 : err += readVal( inrec, String("jwidth"), jwidth );
4076 0 : err += readVal( inrec, String("minweight"), minWeight );
4077 0 : 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 0 : if (ftmachine=="awprojectft" && cfCache=="") {
4083 0 : cfCache = imageName + ".cf";
4084 : }
4085 :
4086 0 : if ( ftmachine=="awprojectft" &&
4087 0 : usePointing==True &&
4088 0 : pointingOffsetSigDev.nelements() != 2 ) {
4089 : // Set the default to a large value so that it behaves like CASA 5.6's usepointing=True.
4090 0 : pointingOffsetSigDev.resize(2);
4091 0 : pointingOffsetSigDev[0] = 600.0;
4092 0 : pointingOffsetSigDev[1] = 600.0;
4093 : }
4094 :
4095 0 : err += verify();
4096 :
4097 0 : } catch(AipsError &x) {
4098 0 : err = err + x.getMesg() + "\n";
4099 0 : }
4100 :
4101 0 : if (err.length()>0) {
4102 0 : throw(AipsError("Invalid Gridding/FTM Parameter set: " + err));
4103 : }
4104 :
4105 0 : }
4106 :
4107 0 : String SynthesisParamsGrid::verify() const
4108 : {
4109 0 : String err;
4110 :
4111 : // Check for valid FTMachine type.
4112 : // Valid other params per FTM type, etc... ( check about nterms>1 )
4113 :
4114 :
4115 0 : if ( imageName == "" ) {
4116 0 : err += "Please supply an image name\n";
4117 : }
4118 0 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
4119 0 : (ftmachine != "mosaicft") && (ftmachine.at(0,3) != "awp") &&
4120 0 : (ftmachine != "mawprojectft") && (ftmachine != "protoft") &&
4121 0 : (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 0 : if ( ( ftmachine == "mosaicft" and mType == "imagemosaic" ) or
4132 0 : ( 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : if ( ftmachine == "awprojectft" and usePointing == True and
4190 0 : 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 0 : if ( ftmachine == "sd" ) {
4198 0 : if ( convertFirst != "always" and
4199 0 : 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 0 : return err;
4207 0 : }
4208 :
4209 0 : void SynthesisParamsGrid::setDefaults()
4210 : {
4211 0 : imageName="";
4212 : // FTMachine parameters
4213 : //ftmachine="GridFT";
4214 0 : ftmachine="gridft";
4215 0 : gridder=ftmachine;
4216 0 : padding=1.2;
4217 0 : useAutoCorr=false;
4218 0 : useDoublePrec=true;
4219 0 : wprojplanes=1;
4220 0 : convFunc="SF";
4221 0 : vpTable="";
4222 :
4223 : // facets
4224 0 : facets=1;
4225 :
4226 : // chanchunks
4227 0 : chanchunks=1;
4228 :
4229 : // Spectral Axis interpolation
4230 0 : interpolation=String("nearest");
4231 :
4232 : //mosaic use pointing
4233 0 : usePointing=false;
4234 : // Moving phase center ?
4235 0 : distance=Quantity(0,"m");
4236 0 : trackSource=false;
4237 0 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
4238 :
4239 : // The extra params for WB-AWP
4240 0 : aTermOn = true;
4241 0 : psTermOn = true;
4242 0 : mTermOn = false;
4243 0 : wbAWP = true;
4244 0 : cfCache = "";
4245 0 : usePointing = false;
4246 0 : pointingOffsetSigDev.resize(0);
4247 : // pointingOffsetSigDev.set(30.0);
4248 0 : doPBCorr = true;
4249 0 : conjBeams = true;
4250 0 : computePAStep=360.0;
4251 0 : rotatePAStep=5.0;
4252 :
4253 : // extra params for single-dish
4254 0 : pointingDirCol = "";
4255 0 : convertFirst = "never";
4256 0 : skyPosThreshold = 0.0;
4257 0 : convSupport = -1;
4258 0 : truncateSize = Quantity(-1.0);
4259 0 : gwidth = Quantity(-1.0);
4260 0 : jwidth = Quantity(-1.0);
4261 0 : minWeight = 0.0;
4262 0 : clipMinMax = False;
4263 :
4264 : // Mapper type
4265 0 : mType = String("default");
4266 :
4267 0 : }
4268 :
4269 0 : Record SynthesisParamsGrid::toRecord() const
4270 : {
4271 0 : Record gridpar;
4272 :
4273 0 : gridpar.define("imagename", imageName);
4274 : // FTMachine params
4275 0 : gridpar.define("padding", padding);
4276 0 : gridpar.define("useautocorr",useAutoCorr );
4277 0 : gridpar.define("usedoubleprec", useDoublePrec);
4278 0 : gridpar.define("wprojplanes", wprojplanes);
4279 0 : gridpar.define("convfunc", convFunc);
4280 0 : gridpar.define("vptable", vpTable);
4281 :
4282 0 : gridpar.define("facets", facets);
4283 0 : gridpar.define("chanchunks", chanchunks);
4284 :
4285 0 : gridpar.define("interpolation",interpolation);
4286 :
4287 0 : gridpar.define("distance", QuantityToString(distance));
4288 0 : gridpar.define("tracksource", trackSource);
4289 0 : gridpar.define("trackdir", MDirectionToString( trackDir ));
4290 :
4291 0 : gridpar.define("aterm",aTermOn );
4292 0 : gridpar.define("psterm",psTermOn );
4293 0 : gridpar.define("mterm",mTermOn );
4294 0 : gridpar.define("wbawp", wbAWP);
4295 0 : gridpar.define("cfcache", cfCache);
4296 0 : gridpar.define("usepointing",usePointing );
4297 0 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
4298 0 : gridpar.define("dopbcorr", doPBCorr);
4299 0 : gridpar.define("conjbeams",conjBeams );
4300 0 : gridpar.define("computepastep", computePAStep);
4301 0 : gridpar.define("rotatepastep", rotatePAStep);
4302 :
4303 0 : gridpar.define("pointingcolumntouse", pointingDirCol );
4304 0 : gridpar.define("convertfirst", convertFirst );
4305 0 : gridpar.define("skyposthreshold", skyPosThreshold );
4306 0 : gridpar.define("convsupport", convSupport );
4307 0 : gridpar.define("truncate", QuantityToString(truncateSize) );
4308 0 : gridpar.define("gwidth", QuantityToString(gwidth) );
4309 0 : gridpar.define("jwidth", QuantityToString(jwidth) );
4310 0 : gridpar.define("minweight", minWeight );
4311 0 : gridpar.define("clipminmax", clipMinMax );
4312 :
4313 0 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
4314 : /// else gridpar.define("deconvolver","singleterm");
4315 :
4316 0 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
4317 0 : else gridpar.define("gridder", gridder);
4318 :
4319 : // gridpar.define("mtype", mType);
4320 :
4321 0 : return gridpar;
4322 0 : }
4323 :
4324 :
4325 :
4326 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
4327 :
4328 : /////////////////////// Deconvolver Parameters
4329 :
4330 0 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
4331 : {
4332 0 : setDefaults();
4333 0 : }
4334 :
4335 0 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
4336 : {
4337 0 : }
4338 :
4339 :
4340 0 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
4341 : {
4342 0 : setDefaults();
4343 :
4344 0 : String err("");
4345 :
4346 : try
4347 : {
4348 :
4349 0 : err += readVal( inrec, String("imagename"), imageName );
4350 0 : 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 0 : if( inrec.isDefined("startmodel") )
4356 : {
4357 0 : if( inrec.dataType("startmodel")==TpString )
4358 : {
4359 0 : String onemodel;
4360 0 : err += readVal( inrec, String("startmodel"), onemodel );
4361 0 : if( onemodel.length()>0 )
4362 : {
4363 0 : startModel.resize(1);
4364 0 : startModel[0] = onemodel;
4365 : }
4366 0 : else {startModel.resize();}
4367 0 : }
4368 0 : else if( inrec.dataType("startmodel")==TpArrayString ||
4369 0 : inrec.dataType("startmodel")==TpArrayBool)
4370 : {
4371 0 : 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 0 : err += readVal( inrec, String("id"), deconvolverId );
4380 0 : err += readVal( inrec, String("nterms"), nTaylorTerms );
4381 :
4382 0 : err += readVal( inrec, String("scales"), scales );
4383 0 : err += readVal( inrec, String("scalebias"), scalebias );
4384 :
4385 0 : err += readVal( inrec, String("usemask"), maskType );
4386 0 : if( maskType=="auto-thresh" )
4387 : {
4388 0 : autoMaskAlgorithm = "thresh";
4389 : }
4390 0 : else if( maskType=="auto-thesh2" )
4391 : {
4392 0 : autoMaskAlgorithm = "thresh2";
4393 : }
4394 0 : else if( maskType=="auto-onebox" )
4395 : {
4396 0 : autoMaskAlgorithm = "onebox";
4397 : }
4398 0 : else if( maskType=="user" || maskType=="pb" )
4399 : {
4400 0 : autoMaskAlgorithm = "";
4401 : }
4402 :
4403 :
4404 0 : if( inrec.isDefined("mask") )
4405 : {
4406 0 : if( inrec.dataType("mask")==TpString )
4407 : {
4408 0 : err+= readVal( inrec, String("mask"), maskString );
4409 : }
4410 0 : else if( inrec.dataType("mask")==TpArrayString )
4411 : {
4412 0 : err+= readVal( inrec, String("mask"), maskList );
4413 : }
4414 : }
4415 :
4416 0 : if( inrec.isDefined("pbmask") )
4417 : {
4418 0 : err += readVal( inrec, String("pbmask"), pbMask );
4419 : }
4420 0 : if( inrec.isDefined("maskthreshold") )
4421 : {
4422 0 : if( inrec.dataType("maskthreshold")==TpString )
4423 : {
4424 0 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
4425 : //deal with the case a string is a float value without unit
4426 0 : Quantity testThresholdString;
4427 0 : Quantity::read(testThresholdString,maskThreshold);
4428 0 : if( testThresholdString.getUnit()=="" )
4429 : {
4430 0 : if(testThresholdString.getValue()<1.0)
4431 : {
4432 0 : fracOfPeak = testThresholdString.getValue();
4433 0 : maskThreshold=String("");
4434 : }
4435 : }
4436 0 : }
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 0 : if( inrec.isDefined("maskresolution") )
4455 : {
4456 0 : if( inrec.dataType("maskresolution")==TpString )
4457 : {
4458 0 : err += readVal(inrec, String("maskresolution"), maskResolution );
4459 : //deal with the case a string is a float value without unit
4460 0 : Quantity testResolutionString;
4461 0 : Quantity::read(testResolutionString,maskResolution);
4462 0 : if( testResolutionString.getUnit()=="" )
4463 : {
4464 0 : maskResByBeam = testResolutionString.getValue();
4465 0 : maskResolution=String("");
4466 : }
4467 0 : }
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 0 : if( inrec.isDefined("nmask") )
4481 : {
4482 0 : if( inrec.dataType("nmask")==TpInt )
4483 : {
4484 0 : err+= readVal(inrec, String("nmask"), nMask );
4485 : }
4486 : else
4487 : {
4488 0 : err+= "nmask must be an integer\n";
4489 : }
4490 : }
4491 0 : if( inrec.isDefined("autoadjust") )
4492 : {
4493 0 : if( inrec.dataType("autoadjust")==TpBool )
4494 : {
4495 0 : 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 0 : if (inrec.isDefined("fusedthreshold"))
4504 : {
4505 0 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
4506 0 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
4507 : else
4508 0 : err += "fusedthreshold must be a float or double";
4509 : }
4510 0 : if (inrec.isDefined("specmode"))
4511 : {
4512 0 : if(inrec.dataType("specmode") == TpString)
4513 0 : 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 0 : if (inrec.isDefined("largestscale"))
4519 : {
4520 0 : if (inrec.dataType("largestscale") == TpInt)
4521 0 : 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 0 : if( inrec.isDefined("sidelobethreshold"))
4527 : {
4528 0 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
4529 : {
4530 0 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
4531 : }
4532 : else
4533 : {
4534 0 : err+= "sidelobethreshold must be a float or double";
4535 : }
4536 : }
4537 :
4538 0 : if( inrec.isDefined("noisethreshold"))
4539 : {
4540 0 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
4541 : {
4542 0 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4543 : }
4544 : else
4545 : {
4546 0 : err+= "noisethreshold must be a float or double";
4547 : }
4548 : }
4549 0 : if( inrec.isDefined("lownoisethreshold"))
4550 : {
4551 0 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4552 : {
4553 0 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4554 : }
4555 : else
4556 : {
4557 0 : err+= "lownoisethreshold must be a float or double";
4558 : }
4559 : }
4560 0 : if( inrec.isDefined("negativethreshold"))
4561 : {
4562 0 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4563 : {
4564 0 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4565 : }
4566 : else
4567 : {
4568 0 : err+= "negativethreshold must be a float or double";
4569 : }
4570 : }
4571 0 : if( inrec.isDefined("smoothfactor"))
4572 : {
4573 0 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4574 : {
4575 0 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4576 : }
4577 : else
4578 : {
4579 0 : err+= "smoothfactor must be a float or double";
4580 : }
4581 : }
4582 0 : if( inrec.isDefined("minbeamfrac"))
4583 : {
4584 0 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4585 : {
4586 0 : 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 0 : if( inrec.isDefined("cutthreshold"))
4600 : {
4601 0 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4602 : {
4603 0 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4604 : }
4605 : else {
4606 0 : err+= "cutthreshold must be a float or double";
4607 : }
4608 : }
4609 0 : if( inrec.isDefined("growiterations"))
4610 : {
4611 0 : if (inrec.dataType("growiterations")==TpInt) {
4612 0 : err+= readVal(inrec, String("growiterations"), growIterations );
4613 : }
4614 : else {
4615 0 : err+= "growiterations must be an integer\n";
4616 : }
4617 : }
4618 0 : if( inrec.isDefined("dogrowprune"))
4619 : {
4620 0 : if (inrec.dataType("dogrowprune")==TpBool) {
4621 0 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4622 : }
4623 : else {
4624 0 : err+= "dogrowprune must be a bool\n";
4625 : }
4626 : }
4627 0 : if( inrec.isDefined("minpercentchange"))
4628 : {
4629 0 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4630 0 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4631 : }
4632 : else {
4633 0 : err+= "minpercentchange must be a float or double";
4634 : }
4635 : }
4636 0 : if( inrec.isDefined("verbose"))
4637 : {
4638 0 : if (inrec.dataType("verbose")==TpBool ) {
4639 0 : err+= readVal(inrec, String("verbose"), verbose);
4640 : }
4641 : else {
4642 0 : err+= "verbose must be a bool";
4643 : }
4644 : }
4645 0 : if( inrec.isDefined("fastnoise"))
4646 : {
4647 0 : if (inrec.dataType("fastnoise")==TpBool ) {
4648 0 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4649 : }
4650 : else {
4651 0 : err+= "fastnoise must be a bool";
4652 : }
4653 : }
4654 0 : if( inrec.isDefined("nsigma") )
4655 : {
4656 0 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4657 0 : 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 0 : if( inrec.isDefined("noRequireSumwt") )
4670 : {
4671 0 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4672 0 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4673 : }
4674 : else {
4675 0 : err+= "noRequireSumwt must be a bool";
4676 : }
4677 : }
4678 0 : if( inrec.isDefined("fullsummary") )
4679 : {
4680 0 : if (inrec.dataType("fullsummary")==TpBool) {
4681 0 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4682 : }
4683 : else {
4684 0 : err+= "fullsummary must be a bool";
4685 : }
4686 : }
4687 0 : if( inrec.isDefined("restoringbeam") )
4688 : {
4689 0 : String errinfo("");
4690 : try {
4691 :
4692 0 : if( inrec.dataType("restoringbeam")==TpString )
4693 : {
4694 0 : 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 0 : if( (! usebeam.matches("common")) && usebeam.length()!=0 )
4699 : {
4700 0 : Quantity bsize;
4701 0 : err += readVal( inrec, String("restoringbeam"), bsize );
4702 0 : restoringbeam.setMajorMinor( bsize, bsize );
4703 0 : usebeam = String("");
4704 0 : }
4705 0 : errinfo = usebeam;
4706 : }
4707 0 : 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 0 : }// if isdefined(restoringbeam)
4741 :
4742 0 : if( inrec.isDefined("interactive") )
4743 : {
4744 0 : if( inrec.dataType("interactive")==TpBool )
4745 0 : {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 0 : err += verify();
4753 :
4754 : }
4755 0 : catch(AipsError &x)
4756 : {
4757 0 : err = err + x.getMesg() + "\n";
4758 0 : }
4759 :
4760 0 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4761 :
4762 0 : }
4763 :
4764 0 : String SynthesisParamsDeconv::verify() const
4765 : {
4766 0 : String err;
4767 :
4768 0 : 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 0 : if( maskString.length()>0 )
4772 : {
4773 0 : File fp( imageName+".mask" );
4774 0 : if( fp.exists() ) err += "Mask image " + imageName+".mask exists, but a specific input mask of " + maskString + " has also been supplied. Please either reset mask='' to reuse the existing mask, or delete " + imageName + ".mask before restarting";
4775 0 : }
4776 :
4777 0 : if( pbMask >= 1.0)
4778 0 : {err += "pbmask must be < 1.0 \n"; }
4779 0 : else if( pbMask < 0.0)
4780 0 : {err += "pbmask must be a positive value \n"; }
4781 :
4782 0 : 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 0 : if ( fracOfPeak >= 1.0)
4791 0 : {err += "fracofpeak must be < 1.0 \n"; }
4792 0 : else if ( fracOfPeak < 0.0)
4793 0 : {err += "fracofpeak must be a positive value \n"; }
4794 :
4795 0 : return err;
4796 0 : }
4797 :
4798 0 : void SynthesisParamsDeconv::setDefaults()
4799 : {
4800 0 : imageName="";
4801 0 : algorithm="hogbom";
4802 0 : startModel=Vector<String>(0);
4803 0 : deconvolverId=0;
4804 0 : nTaylorTerms=1;
4805 0 : scales.resize(1); scales[0]=0.0;
4806 0 : scalebias=0.6;
4807 0 : maskType="none";
4808 0 : maskString="";
4809 0 : maskList.resize(1); maskList[0]="";
4810 0 : pbMask=0.0;
4811 0 : autoMaskAlgorithm="thresh";
4812 0 : maskThreshold="";
4813 0 : maskResolution="";
4814 0 : fracOfPeak=0.0;
4815 0 : nMask=0;
4816 0 : interactive=false;
4817 0 : autoAdjust=False;
4818 0 : fusedThreshold = 0.0;
4819 0 : specmode="mfs";
4820 0 : largestscale = -1;
4821 0 : }
4822 :
4823 0 : Record SynthesisParamsDeconv::toRecord() const
4824 : {
4825 0 : Record decpar;
4826 :
4827 0 : decpar.define("imagename", imageName);
4828 0 : decpar.define("deconvolver", algorithm);
4829 0 : decpar.define("startmodel",startModel);
4830 0 : decpar.define("id",deconvolverId);
4831 0 : decpar.define("nterms",nTaylorTerms);
4832 0 : decpar.define("scales",scales);
4833 0 : decpar.define("scalebias",scalebias);
4834 0 : decpar.define("usemask",maskType);
4835 0 : decpar.define("fusedthreshold", fusedThreshold);
4836 0 : decpar.define("specmode", specmode);
4837 0 : decpar.define("largestscale", largestscale);
4838 0 : if( maskList.nelements()==1 && maskList[0]=="")
4839 : {
4840 0 : decpar.define("mask",maskString);
4841 : }
4842 : else {
4843 0 : decpar.define("mask",maskList);
4844 : }
4845 0 : decpar.define("pbmask",pbMask);
4846 0 : if (fracOfPeak > 0.0)
4847 : {
4848 0 : decpar.define("maskthreshold",fracOfPeak);
4849 : }
4850 : else
4851 : {
4852 0 : decpar.define("maskthreshold",maskThreshold);
4853 : }
4854 0 : decpar.define("maskresolution",maskResolution);
4855 0 : decpar.define("nmask",nMask);
4856 0 : decpar.define("autoadjust",autoAdjust);
4857 0 : decpar.define("sidelobethreshold",sidelobeThreshold);
4858 0 : decpar.define("noisethreshold",noiseThreshold);
4859 0 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4860 0 : decpar.define("negativethreshold",negativeThreshold);
4861 0 : decpar.define("smoothfactor",smoothFactor);
4862 0 : decpar.define("minbeamfrac",minBeamFrac);
4863 0 : decpar.define("cutthreshold",cutThreshold);
4864 0 : decpar.define("growiterations",growIterations);
4865 0 : decpar.define("dogrowprune",doGrowPrune);
4866 0 : decpar.define("minpercentchange",minPercentChange);
4867 0 : decpar.define("verbose", verbose);
4868 0 : decpar.define("fastnoise", fastnoise);
4869 0 : decpar.define("interactive",interactive);
4870 0 : decpar.define("nsigma",nsigma);
4871 0 : decpar.define("noRequireSumwt",noRequireSumwt);
4872 0 : decpar.define("fullsummary",fullsummary);
4873 :
4874 0 : return decpar;
4875 0 : }
4876 :
4877 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4878 :
4879 :
4880 : } //# NAMESPACE CASA - END
4881 :
|
