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