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