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