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 0 : SynthesisUtilMethods::SynthesisUtilMethods()
86 : {
87 :
88 0 : }
89 :
90 0 : SynthesisUtilMethods::~SynthesisUtilMethods()
91 : {
92 0 : }
93 :
94 0 : Int SynthesisUtilMethods::validate(const VisBuffer& vb)
95 : {
96 0 : Int N=vb.nRow(),M=-1;
97 0 : for(Int i=0;i<N;i++)
98 : {
99 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
100 0 : {M++;break;}
101 : }
102 0 : return M;
103 : }
104 :
105 0 : Int SynthesisUtilMethods::validate(const vi::VisBuffer2& vb)
106 : {
107 0 : Int N=vb.nRows(),M=-1;
108 0 : for(Int i=0;i<N;i++)
109 : {
110 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
111 0 : {M++;break;}
112 : }
113 0 : return M;
114 : }
115 : // Get the next largest even composite of 2,3,5.
116 : // This is to ensure a 'good' image size for FFTW.
117 : // Translated from gcwrap/scripts/cleanhelper.py : getOptimumSize
118 0 : Int SynthesisUtilMethods::getOptimumSize(const Int npix)
119 : {
120 0 : Int n=npix;
121 :
122 0 : if( n%2 !=0 ){ n+= 1; }
123 :
124 0 : Vector<uInt> fac = primeFactors(n, false);
125 : Int val, newlarge;
126 0 : for( uInt k=0; k< fac.nelements(); k++ )
127 : {
128 0 : if( fac[k]>5 )
129 : {
130 0 : val = fac[k];
131 0 : while( max( primeFactors(val) ) > 5 ){ val+=1;}
132 0 : fac[k] = val;
133 : }
134 : }
135 0 : newlarge=product(fac);
136 0 : for( Int k=n; k<newlarge; k+=2 )
137 : {
138 0 : if( max( primeFactors(k) ) < 6 ) {return k;}
139 : }
140 0 : return newlarge;
141 : }
142 :
143 : // Return the prime factors of the given number
144 0 : Vector<uInt> SynthesisUtilMethods::primeFactors(uInt n, Bool /*douniq*/)
145 : {
146 0 : Vector<uInt> factors;
147 :
148 0 : Int lastresult = n;
149 0 : Int sqlast = int(sqrt(n))+1;
150 :
151 0 : if(n==1){ factors.resize(1);factors[0]=1;return factors;}
152 0 : Int c=2;
153 : while(1)
154 : {
155 0 : if( lastresult == 1 || c > sqlast ) { break; }
156 0 : sqlast = (Int)(sqrt(lastresult))+1;
157 : while(1)
158 : {
159 0 : if( c>sqlast){ c=lastresult; break; }
160 0 : if( lastresult % c == 0 ) { break; }
161 0 : c += 1;
162 : }
163 0 : factors.resize( factors.nelements()+1, true );
164 0 : factors[ factors.nelements()-1 ] = c;
165 0 : lastresult /= c;
166 : }
167 0 : if( factors.nelements()==0 ) { factors.resize(1);factors[0]=n; }
168 :
169 : //if( douniq ) { factors = unique(factors); }
170 :
171 : /*
172 : /// The Sort::unique isn't working as called below. Results appear to be the
173 : /// same as with the cleanhelper python code, so leaving as is for not. CAS-7889
174 : if( douniq )
175 : {
176 : cout << "Test unique fn on : " << factors << endl;
177 : Sort srt;
178 : Vector<uInt> unvec=factors; uInt nrec;
179 : srt.unique(unvec,nrec);
180 : cout << " Unique : " << unvec << " nr : " << nrec << endl;
181 : }
182 : */
183 :
184 0 : return factors;
185 : }
186 :
187 :
188 0 : Bool SynthesisUtilMethods::fitPsfBeam(const String& imagename, const Int nterms, const Float psfcutoff)
189 : {
190 0 : LogIO os(LogOrigin("SynthesisUtilMethods", "fitPsfBeam"));
191 :
192 0 : if (psfcutoff >=1.0 || psfcutoff<=0.0)
193 : {
194 0 : os << "psfcutoff must be >0 and <1" << LogIO::WARN;
195 0 : return false;
196 : }
197 :
198 0 : std::shared_ptr<SIImageStore> imstore;
199 0 : if( nterms>1 )
200 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStoreMultiTerm( imagename, nterms, true, true )); }
201 : else
202 0 : { imstore = std::shared_ptr<SIImageStore>(new SIImageStore( imagename, true, true )); }
203 :
204 :
205 0 : os << "Fitting PSF beam for Imagestore : " << imstore->getName() << LogIO::POST;
206 :
207 0 : imstore->makeImageBeamSet(psfcutoff, true);
208 :
209 0 : imstore->printBeamSet();
210 :
211 0 : imstore->releaseLocks();
212 :
213 0 : return true;
214 : }
215 :
216 :
217 :
218 :
219 :
220 : /***make a record of synthesisimager::weight parameters***/
221 0 : 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 0 : Record outRec;
229 0 : outRec.define("type", type);
230 0 : outRec.define("rmode", rmode);
231 0 : Record quantRec;
232 0 : QuantumHolder(noise).toRecord(quantRec);
233 0 : outRec.defineRecord("noise", quantRec);
234 0 : outRec.define("robust", robust);
235 0 : QuantumHolder(fieldofview).toRecord(quantRec);
236 0 : outRec.defineRecord("fieldofview", quantRec);
237 0 : outRec.define("npixels", npixels);
238 0 : outRec.define("multifield", multiField);
239 0 : outRec.define("usecubebriggs", useCubeBriggs);
240 0 : outRec.define("filtertype", filtertype);
241 0 : QuantumHolder(filterbmaj).toRecord(quantRec);
242 0 : outRec.defineRecord("filterbmaj", quantRec);
243 0 : QuantumHolder(filterbmin).toRecord(quantRec);
244 0 : outRec.defineRecord("filterbmin", quantRec);
245 0 : QuantumHolder(filterbpa).toRecord(quantRec);
246 0 : outRec.defineRecord("filterbpa", quantRec);
247 0 : outRec.define("fracBW", fracBW);
248 :
249 0 : return outRec;
250 : }
251 0 : 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 0 : QuantumHolder qh;
258 0 : String err;
259 0 : if(!inRec.isDefined("type"))
260 0 : throw(AipsError("Record is not filled with SynthesisUtilMethods::fillWeightRecord"));
261 0 : inRec.get("type", type);
262 0 : inRec.get("rmode", rmode);
263 0 : if(!qh.fromRecord(err, inRec.asRecord("noise")))
264 0 : throw(AipsError("Error in reading noise param"));
265 0 : noise=qh.asQuantity();
266 0 : inRec.get("robust", robust);
267 0 : if(!qh.fromRecord(err, inRec.asRecord("fieldofview")))
268 0 : throw(AipsError("Error in reading fieldofview param"));
269 0 : fieldofview=qh.asQuantity();
270 0 : inRec.get("npixels", npixels);
271 0 : inRec.get("multifield", multiField);
272 0 : inRec.get("usecubebriggs", useCubeBriggs);
273 0 : inRec.get("filtertype", filtertype);
274 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmaj")))
275 0 : throw(AipsError("Error in reading filterbmaj param"));
276 0 : filterbmaj=qh.asQuantity();
277 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbmin")))
278 0 : throw(AipsError("Error in reading filterbmin param"));
279 0 : filterbmin=qh.asQuantity();
280 0 : if(!qh.fromRecord(err, inRec.asRecord("filterbpa")))
281 0 : throw(AipsError("Error in reading filterbpa param"));
282 0 : filterbpa=qh.asQuantity();
283 0 : inRec.get("fracBW", fracBW);
284 :
285 :
286 :
287 0 : }
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 0 : Int SynthesisUtilMethods::parseProcStatusLine(const std::string &str) {
297 0 : istringstream is(str);
298 0 : std::string token;
299 0 : is >> token;
300 0 : is >> token;
301 0 : Int value = stoi(token);
302 0 : 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 0 : String SynthesisUtilMethods::makeResourceFilename(int pid)
315 : {
316 0 : 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 0 : const int strMax = 255;
320 : char hostname[strMax];
321 0 : gethostname(hostname, strMax);
322 0 : g_hostname = hostname;
323 :
324 0 : auto time = std::time(nullptr);
325 0 : auto gmt = std::gmtime(&time);
326 0 : const char* format = "%Y%m%d_%H%M%S";
327 : char timestr[strMax];
328 0 : std::strftime(timestr, strMax, format, gmt);
329 0 : g_startTimestamp = timestr;
330 : }
331 :
332 0 : return String("casa.synthesis.imager.memprofile." + String::toString(pid) +
333 0 : "." + 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 0 : 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 0 : Bool isOn = false;
448 0 : AipsrcValue<Bool>::find(isOn, g_enableOptMemProfile);
449 0 : if (!isOn)
450 0 : return;
451 :
452 : // TODO: reorganize, use struct or something to hold and pass info over. ifdef lnx
453 0 : LogIO casalog( LogOrigin("SynthesisUtilMethods", "getResource", WHERE) );
454 :
455 : // To hold memory stats, in MB
456 0 : int vmRSS = -1, vmHWM = -1, vmSize = -1, vmPeak = -1, vmSwap = -1;
457 0 : pid_t pid = -1;
458 0 : int fdSize = -1;
459 :
460 : // TODO: this won't probably work on anything but linux
461 0 : ifstream procFile("/proc/self/status");
462 0 : if (procFile.is_open()) {
463 0 : std::string line;
464 0 : while (not procFile.eof()) {
465 0 : getline(procFile, line);
466 0 : const std::string startVmRSS = "VmRSS:";
467 0 : const std::string startVmWHM = "VmHWM:";
468 0 : const std::string startVmSize = "VmSize:";
469 0 : const std::string startVmPeak = "VmPeak:";
470 0 : const std::string startVmSwap = "VmSwap:";
471 0 : const std::string startFDSize = "FDSize:";
472 0 : const double KB_TO_MB = 1024.0;
473 0 : if (startVmRSS == line.substr(0, startVmRSS.size())) {
474 0 : vmRSS = parseProcStatusLine(line.c_str()) / KB_TO_MB;
475 0 : } else if (startVmWHM == line.substr(0, startVmWHM.size())) {
476 0 : vmHWM = parseProcStatusLine(line.c_str()) / KB_TO_MB;
477 0 : } else if (startVmSize == line.substr(0, startVmSize.size())) {
478 0 : vmSize = parseProcStatusLine(line.c_str()) / KB_TO_MB;
479 0 : } else if (startVmPeak == line.substr(0, startVmPeak.size())) {
480 0 : vmPeak = parseProcStatusLine(line.c_str()) / KB_TO_MB;
481 0 : } else if (startVmSwap == line.substr(0, startVmSwap.size())) {
482 0 : vmSwap = parseProcStatusLine(line.c_str()) / KB_TO_MB;
483 0 : } else if (startFDSize == line.substr(0, startFDSize.size())) {
484 0 : fdSize = parseProcStatusLine(line.c_str());
485 : }
486 : }
487 0 : procFile.close();
488 : }
489 :
490 0 : pid = getpid();
491 :
492 : struct rusage usage;
493 : struct timeval now;
494 0 : getrusage(RUSAGE_SELF, &usage);
495 0 : 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 0 : if (vmHWM < 0) {
501 0 : vmHWM = usage.ru_maxrss;
502 : }
503 :
504 0 : ostringstream oss;
505 0 : oss << "PID: " << pid ;
506 0 : oss << " MemRSS (VmRSS): " << vmRSS << " MB.";
507 0 : oss << " VmWHM: " << vmHWM << " MB.";
508 0 : oss << " VirtMem (VmSize): " << vmSize << " MB.";
509 0 : oss << " VmPeak: " << vmPeak << " MB.";
510 0 : oss << " VmSwap: " << vmSwap << " MB.";
511 0 : oss << " ProcTime: " << now.tv_sec << '.' << now.tv_usec;
512 0 : oss << " FDSize: " << fdSize;
513 0 : oss << " [" << label << "] ";
514 0 : casalog << oss.str() << LogIO::NORMAL3 << LogIO::POST;
515 :
516 : // Write this to a file too...
517 : try {
518 0 : if (fname.empty()) {
519 0 : fname = makeResourceFilename(pid);
520 : }
521 0 : ofstream ofile(fname, ios::app);
522 0 : if (ofile.is_open()) {
523 0 : if (0 == ofile.tellp()) {
524 : casalog << g_enableOptMemProfile << " is enabled, initializing output file for "
525 : "imager profiling information (memory and run time): " << fname <<
526 0 : LogIO::NORMAL << LogIO::POST;
527 0 : ostringstream header;
528 : header << "# PID, MemRSS_(VmRSS)_MB, VmWHM_MB, VirtMem_(VmSize)_MB, VmPeak_MB, "
529 0 : "VmSwap_MB, ProcTime_sec, FDSize, label_checkpoint";
530 0 : ofile << header.str() << '\n';
531 : }
532 0 : ostringstream line;
533 0 : line << pid << ',' << vmRSS << ',' << vmHWM << ',' << vmSize << ','
534 0 : << vmPeak << ','<< vmSwap << ',' << now.tv_sec << '.' << now.tv_usec << ','
535 0 : << fdSize << ',' << '[' << label << ']';
536 0 : ofile << line.str() << '\n';
537 0 : 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 0 : String SynthesisUtilMethods::asComprehensibleDirectionString(MDirection const &direction)
922 : {
923 0 : MVAngle mvRa=direction.getAngle().getValue()(0);
924 0 : MVAngle mvDec=direction.getAngle().getValue()(1);
925 0 : ostringstream oos;
926 0 : oos << " ";
927 0 : Int widthRA=20;
928 0 : Int widthDec=20;
929 0 : oos.setf(ios::left, ios::adjustfield);
930 0 : oos.width(widthRA); oos << mvRa(0.0).string(MVAngle::TIME,8);
931 0 : oos.width(widthDec); oos << mvDec.string(MVAngle::DIG2,8);
932 0 : oos << " "
933 0 : << MDirection::showType(direction.getRefPtr()->getType());
934 0 : return String(oos);
935 : }
936 :
937 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
938 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
939 : ///////////////// Parameter Containers ///////////////////////////////////////////////////////
940 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
941 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
942 :
943 : // Read String from Record
944 0 : String SynthesisParams::readVal(const Record &rec, String id, String& val) const
945 : {
946 0 : if( rec.isDefined( id ) )
947 : {
948 0 : String inval("");
949 0 : if( rec.dataType( id )==TpString )
950 0 : { 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 0 : if(inval.length()>0){val=inval;}
956 0 : return String("");
957 : }
958 0 : else { return String(id + " must be a string\n"); }
959 : }
960 0 : else { return String("");}
961 : }
962 :
963 : // Read Integer from Record
964 0 : String SynthesisParams::readVal(const Record &rec, String id, Int& val) const
965 : {
966 0 : if( rec.isDefined( id ) )
967 : {
968 0 : if( rec.dataType( id )==TpInt ) { rec.get( id , val ); return String(""); }
969 0 : else { return String(id + " must be an integer\n"); }
970 : }
971 0 : else { return String(""); }
972 : }
973 :
974 : // Read Float from Record
975 0 : String SynthesisParams::readVal(const Record &rec, String id, Float& val) const
976 : {
977 0 : if( rec.isDefined( id ) )
978 : {
979 0 : if ( rec.dataType( id )==TpFloat || rec.dataType( id )==TpDouble )
980 0 : { rec.get( id , val ); return String(""); }
981 0 : else { return String(id + " must be a float\n"); }
982 : }
983 0 : else { return String(""); }
984 : }
985 :
986 : // Read Bool from Record
987 0 : String SynthesisParams::readVal(const Record &rec, String id, Bool& val) const
988 : {
989 0 : if( rec.isDefined( id ) )
990 : {
991 0 : if( rec.dataType( id )==TpBool ) { rec.get( id , val ); return String(""); }
992 0 : else { return String(id + " must be a bool\n"); }
993 : }
994 0 : else{ return String(""); }
995 : }
996 :
997 : // Read Vector<Int> from Record
998 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Int>& val) const
999 : {
1000 0 : if( rec.isDefined( id ) )
1001 : {
1002 0 : if( rec.dataType( id )==TpArrayInt || rec.dataType( id )==TpArrayUInt )
1003 0 : { 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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<Float>& val) const
1017 : {
1018 0 : if( rec.isDefined( id ) )
1019 : {
1020 0 : if( rec.dataType( id )==TpArrayFloat )
1021 : {
1022 0 : 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 0 : else if ( rec.dataType( id ) ==TpArrayDouble )
1036 : {
1037 0 : Vector<Double> vec; rec.get( id, vec );
1038 0 : val.resize(vec.nelements());
1039 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1040 0 : return String("");
1041 : }
1042 0 : else if ( rec.dataType( id ) ==TpArrayInt )
1043 : {
1044 0 : Vector<Int> vec; rec.get( id, vec );
1045 0 : val.resize(vec.nelements());
1046 0 : for(uInt i=0;i<val.nelements();i++){val[i]=(Float)vec[i];}
1047 0 : return String("");
1048 : }
1049 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1050 : {
1051 0 : Vector<Bool> vec; rec.get( id, vec );
1052 0 : 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 0 : String SynthesisParams::readVal(const Record &rec, String id, Vector<String>& val) const
1063 : {
1064 0 : if( rec.isDefined( id ) )
1065 : {
1066 0 : if( rec.dataType( id )==TpArrayString || rec.dataType( id )==TpArrayChar )
1067 0 : { rec.get( id , val ); return String(""); }
1068 0 : else if ( rec.dataType( id ) == TpArrayBool ) // An empty python vector [] comes in as this.
1069 : {
1070 0 : Vector<Bool> vec; rec.get( id, vec );
1071 0 : 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 0 : 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 0 : 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 0 : String SynthesisParams::stringToMDirection(String instr, MDirection& md) const
1094 : {
1095 : try
1096 : {
1097 0 : String tmpRF, tmpRA, tmpDEC;
1098 0 : std::istringstream iss(instr);
1099 0 : iss >> tmpRF >> tmpRA >> tmpDEC;
1100 0 : 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 0 : casacore::Quantity tmpQRA;
1107 0 : casacore::Quantity tmpQDEC;
1108 0 : casacore::Quantity::read(tmpQRA, tmpRA);
1109 0 : casacore::Quantity::read(tmpQDEC, tmpDEC);
1110 :
1111 0 : 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 0 : Bool status = MDirection::getType(theRF, tmpRF);
1122 0 : if (!status) {
1123 0 : throw AipsError();
1124 : }
1125 0 : md = MDirection (tmpQRA, tmpQDEC, theRF);
1126 0 : return String("");
1127 : }
1128 0 : catch(AipsError &x)
1129 : {
1130 0 : 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 0 : String SynthesisParams::readVal(const Record &rec, String id, Quantity& val) const
1136 : {
1137 0 : if( rec.isDefined( id ) )
1138 : {
1139 0 : if( rec.dataType( id )==TpString )
1140 0 : { 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 0 : else{ return String(""); }
1144 : }
1145 :
1146 : // Read MDirection from a Record string
1147 0 : String SynthesisParams::readVal(const Record &rec, String id, MDirection& val) const
1148 : {
1149 0 : if( rec.isDefined( id ) )
1150 : {
1151 0 : if( rec.dataType( id )==TpString )
1152 0 : { 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 0 : else{ return String(""); }
1156 : }
1157 :
1158 : // Convert MDirection to String
1159 0 : String SynthesisParams::MDirectionToString(MDirection val) const
1160 : {
1161 0 : MVDirection mvpc( val.getAngle() );
1162 0 : MVAngle ra = mvpc.get()(0);
1163 0 : MVAngle dec = mvpc.get()(1);
1164 : // Beware of precision here......( for VLBA / ALMA ). 14 gives 8 digits after the decimal for arcsec.
1165 0 : return String(val.getRefString() + " " + ra.string(MVAngle::TIME,14) + " " + dec.string(MVAngle::ANGLE,14));
1166 : }
1167 :
1168 : // Convert Quantity to String
1169 0 : String SynthesisParams::QuantityToString(Quantity val) const
1170 : {
1171 0 : 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 0 : ss.precision(std::numeric_limits<double>::digits10+2);
1176 0 : ss << val;
1177 0 : 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 0 : String SynthesisParams::recordQMToString(const Record &rec) const
1186 : {
1187 : Double val;
1188 0 : String unit;
1189 0 : if ( rec.isDefined("m0") )
1190 : {
1191 0 : Record subrec = rec.subRecord("m0");
1192 0 : subrec.get("value",val);
1193 0 : subrec.get("unit",unit);
1194 : }
1195 0 : else if (rec.isDefined("value") )
1196 : {
1197 0 : rec.get("value",val);
1198 0 : rec.get("unit",unit);
1199 : }
1200 0 : return String::toString(val) + unit;
1201 : }
1202 :
1203 :
1204 : /////////////////////// Selection Parameters
1205 :
1206 0 : SynthesisParamsSelect::SynthesisParamsSelect():SynthesisParams()
1207 : {
1208 0 : setDefaults();
1209 0 : }
1210 :
1211 0 : SynthesisParamsSelect::~SynthesisParamsSelect()
1212 : {
1213 0 : }
1214 :
1215 0 : SynthesisParamsSelect::SynthesisParamsSelect(const SynthesisParamsSelect& other){
1216 0 : operator=(other);
1217 0 : }
1218 :
1219 0 : SynthesisParamsSelect& SynthesisParamsSelect::operator=(const SynthesisParamsSelect& other){
1220 0 : if(this!=&other) {
1221 0 : msname=other.msname;
1222 0 : spw=other.spw;
1223 0 : freqbeg=other.freqbeg;
1224 0 : freqend=other.freqend;
1225 0 : freqframe=other.freqframe;
1226 0 : field=other.field;
1227 0 : antenna=other.antenna;
1228 0 : timestr=other.timestr;
1229 0 : scan=other.scan;
1230 0 : obs=other.obs;
1231 0 : state=other.state;
1232 0 : uvdist=other.uvdist;
1233 0 : taql=other.taql;
1234 0 : usescratch=other.usescratch;
1235 0 : readonly=other.readonly;
1236 0 : incrmodel=other.incrmodel;
1237 0 : datacolumn=other.datacolumn;
1238 :
1239 : }
1240 0 : return *this;
1241 : }
1242 :
1243 0 : void SynthesisParamsSelect::fromRecord(const Record &inrec)
1244 : {
1245 0 : setDefaults();
1246 :
1247 0 : String err("");
1248 :
1249 : try
1250 : {
1251 :
1252 0 : err += readVal( inrec, String("msname"), msname );
1253 :
1254 0 : err += readVal( inrec, String("readonly"), readonly );
1255 0 : err += readVal( inrec, String("usescratch"), usescratch );
1256 :
1257 : // override with entries from savemodel.
1258 0 : String savemodel("");
1259 0 : err += readVal( inrec, String("savemodel"), savemodel );
1260 0 : if( savemodel=="none" ){usescratch=false; readonly=true;}
1261 0 : else if( savemodel=="virtual" ){usescratch=false; readonly=false;}
1262 0 : else if ( savemodel=="modelcolumn" ){ usescratch=true; readonly=false; }
1263 :
1264 0 : err += readVal( inrec, String("incrmodel"), incrmodel );
1265 :
1266 0 : err += readVal( inrec, String("spw"), spw );
1267 :
1268 : /// -------------------------------------------------------------------------------------
1269 : // Why are these params here ? Repeats of defineimage.
1270 0 : err += readVal( inrec, String("freqbeg"), freqbeg);
1271 0 : err += readVal( inrec, String("freqend"), freqend);
1272 :
1273 0 : String freqframestr( MFrequency::showType(freqframe) );
1274 0 : err += readVal( inrec, String("outframe"), freqframestr);
1275 0 : if( ! MFrequency::getType(freqframe, freqframestr) )
1276 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1277 : /// -------------------------------------------------------------------------------------
1278 :
1279 0 : err += readVal( inrec, String("field"),field);
1280 0 : err += readVal( inrec, String("antenna"),antenna);
1281 0 : err += readVal( inrec, String("timestr"),timestr);
1282 0 : err += readVal( inrec, String("scan"),scan);
1283 0 : err += readVal( inrec, String("obs"),obs);
1284 0 : err += readVal( inrec, String("state"),state);
1285 0 : err += readVal( inrec, String("uvdist"),uvdist);
1286 0 : err += readVal( inrec, String("taql"),taql);
1287 :
1288 0 : err += readVal( inrec, String("datacolumn"),datacolumn);
1289 :
1290 0 : err += verify();
1291 :
1292 : }
1293 0 : catch(AipsError &x)
1294 : {
1295 0 : err = err + x.getMesg() + "\n";
1296 : }
1297 :
1298 0 : if( err.length()>0 ) throw(AipsError("Invalid Selection Parameter set : " + err));
1299 :
1300 0 : }
1301 :
1302 0 : String SynthesisParamsSelect::verify() const
1303 : {
1304 0 : String err;
1305 : // Does the MS exist on disk.
1306 0 : Directory thems( msname );
1307 0 : if( thems.exists() )
1308 : {
1309 : // Is it readable ?
1310 0 : if( ! thems.isReadable() )
1311 0 : { err += "MS " + msname + " is not readable.\n"; }
1312 : // Depending on 'readonly', is the MS writable ?
1313 0 : 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 0 : return err;
1320 : }
1321 :
1322 0 : void SynthesisParamsSelect::setDefaults()
1323 : {
1324 0 : msname="";
1325 0 : spw="";
1326 0 : freqbeg="";
1327 0 : freqend="";
1328 0 : MFrequency::getType(freqframe,"LSRK");
1329 0 : field="";
1330 0 : antenna="";
1331 0 : timestr="";
1332 0 : scan="";
1333 0 : obs="";
1334 0 : state="";
1335 0 : uvdist="";
1336 0 : taql="";
1337 0 : usescratch=false;
1338 0 : readonly=true;
1339 0 : incrmodel=false;
1340 0 : datacolumn="corrected";
1341 0 : }
1342 :
1343 0 : Record SynthesisParamsSelect::toRecord()const
1344 : {
1345 0 : Record selpar;
1346 0 : selpar.define("msname",msname);
1347 0 : selpar.define("spw",spw);
1348 0 : selpar.define("freqbeg",freqbeg);
1349 0 : selpar.define("freqend",freqend);
1350 0 : selpar.define("freqframe", MFrequency::showType(freqframe)); // Convert MFrequency::Types to String
1351 : //looks like fromRecord looks for outframe !
1352 0 : selpar.define("outframe", MFrequency::showType(freqframe));
1353 0 : selpar.define("field",field);
1354 0 : selpar.define("antenna",antenna);
1355 0 : selpar.define("timestr",timestr);
1356 0 : selpar.define("scan",scan);
1357 0 : selpar.define("obs",obs);
1358 0 : selpar.define("state",state);
1359 0 : selpar.define("uvdist",uvdist);
1360 0 : selpar.define("taql",taql);
1361 0 : selpar.define("usescratch",usescratch);
1362 0 : selpar.define("readonly",readonly);
1363 0 : selpar.define("incrmodel",incrmodel);
1364 0 : selpar.define("datacolumn",datacolumn);
1365 :
1366 0 : return selpar;
1367 : }
1368 :
1369 :
1370 : /////////////////////// Image Parameters
1371 :
1372 0 : SynthesisParamsImage::SynthesisParamsImage():SynthesisParams()
1373 : {
1374 0 : setDefaults();
1375 0 : }
1376 :
1377 0 : SynthesisParamsImage::~SynthesisParamsImage()
1378 : {
1379 0 : }
1380 :
1381 0 : SynthesisParamsImage& SynthesisParamsImage::operator=(const SynthesisParamsImage& other){
1382 0 : if(this != &other){
1383 0 : imageName=other.imageName;
1384 0 : stokes=other.stokes;
1385 0 : startModel.resize(); startModel=other.startModel;
1386 0 : imsize.resize(); imsize=other.imsize;
1387 0 : cellsize.resize(); cellsize=other.cellsize;
1388 0 : projection=other.projection;
1389 0 : useNCP=other.useNCP;
1390 0 : phaseCenter=other.phaseCenter;
1391 0 : phaseCenterFieldId=other.phaseCenterFieldId;
1392 0 : obslocation=other.obslocation;
1393 0 : pseudoi=other.pseudoi;
1394 0 : nchan=other.nchan;
1395 0 : nTaylorTerms=other.nTaylorTerms;
1396 0 : chanStart=other.chanStart;
1397 0 : chanStep=other.chanStep;
1398 0 : freqStart=other.freqStart;
1399 0 : freqStep=other.freqStep;
1400 0 : refFreq=other.refFreq;
1401 0 : velStart=other.velStart;
1402 0 : velStep=other.velStep;
1403 0 : freqFrame=other.freqFrame;
1404 0 : mFreqStart=other.mFreqStart;
1405 0 : mFreqStep=other.mFreqStep;
1406 0 : mVelStart=other.mVelStart;
1407 0 : mVelStep=other.mVelStep;
1408 0 : restFreq.resize(); restFreq=other.restFreq;
1409 0 : start=other.start;
1410 0 : step=other.step;
1411 0 : frame=other.frame;
1412 0 : veltype=other.veltype;
1413 0 : mode=other.mode;
1414 0 : reffreq=other.reffreq;
1415 0 : sysvel=other.sysvel;
1416 0 : sysvelframe=other.sysvelframe;
1417 0 : sysvelvalue=other.sysvelvalue;
1418 0 : qmframe=other.qmframe;
1419 0 : mveltype=other.mveltype;
1420 0 : tststr=other.tststr;
1421 0 : startRecord=other.startRecord;
1422 0 : stepRecord=other.stepRecord;
1423 0 : reffreqRecord=other.reffreqRecord;
1424 0 : sysvelRecord=other.sysvelRecord;
1425 0 : restfreqRecord=other.restfreqRecord;
1426 0 : csysRecord=other.csysRecord;
1427 0 : csys=other.csys;
1428 0 : imshape.resize(); imshape=other.imshape;
1429 0 : freqFrameValid=other.freqFrameValid;
1430 0 : overwrite=other.overwrite;
1431 0 : deconvolver=other.deconvolver;
1432 0 : distance=other.distance;
1433 0 : trackDir=other.trackDir;
1434 0 : trackSource=other.trackSource;
1435 0 : movingSource=other.movingSource;
1436 :
1437 :
1438 :
1439 : }
1440 :
1441 0 : return *this;
1442 :
1443 : }
1444 :
1445 0 : void SynthesisParamsImage::fromRecord(const Record &inrec)
1446 : {
1447 0 : setDefaults();
1448 0 : String err("");
1449 :
1450 : try
1451 : {
1452 :
1453 0 : err += readVal( inrec, String("imagename"), imageName);
1454 :
1455 : //// imsize
1456 0 : if( inrec.isDefined("imsize") )
1457 : {
1458 0 : DataType tp = inrec.dataType("imsize");
1459 :
1460 0 : if( tp == TpInt ) // A single integer for both dimensions.
1461 : {
1462 0 : Int npix; inrec.get("imsize", npix);
1463 0 : imsize.resize(2);
1464 0 : imsize.set( npix );
1465 : }
1466 0 : else if( tp == TpArrayInt ) // An integer array : [ nx ] or [ nx, ny ]
1467 : {
1468 0 : Vector<Int> ims;
1469 0 : inrec.get("imsize", ims);
1470 0 : if( ims.nelements()==1 ) // [ nx ]
1471 0 : {imsize.set(ims[0]); }
1472 0 : else if( ims.nelements()==2 ) // [ nx, ny ]
1473 0 : { 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 0 : if( inrec.isDefined("cell") )
1484 : {
1485 : try
1486 : {
1487 0 : DataType tp = inrec.dataType("cell");
1488 0 : if( tp == TpInt ||
1489 0 : tp == TpFloat ||
1490 : tp == TpDouble )
1491 : {
1492 0 : Double cell = inrec.asDouble("cell");
1493 0 : cellsize.set( Quantity( cell , "arcsec" ) );
1494 : }
1495 0 : else if ( tp == TpArrayInt ||
1496 0 : tp == TpArrayFloat ||
1497 : tp == TpArrayDouble )
1498 : {
1499 0 : Vector<Double> cells;
1500 0 : inrec.get("cell", cells);
1501 0 : if(cells.nelements()==1) // [ cellx ]
1502 0 : {cellsize.set( Quantity( cells[0], "arcsec" ) ); }
1503 0 : else if( cells.nelements()==2 ) // [ cellx, celly ]
1504 0 : { cellsize[0]=Quantity(cells[0],"arcsec"); cellsize[1]=Quantity(cells[1],"arcsec"); }
1505 : else // Wrong array length
1506 0 : {err += "cellsize must be a single integer/string, or a vector of two integers/strings\n"; }
1507 : }
1508 0 : else if( tp == TpString )
1509 : {
1510 0 : String cell;
1511 0 : inrec.get("cell",cell);
1512 0 : Quantity qcell;
1513 0 : err += stringToQuantity( cell, qcell );
1514 0 : cellsize.set( qcell );
1515 : }
1516 0 : else if( tp == TpArrayString )
1517 : {
1518 0 : Array<String> cells;
1519 0 : inrec.get("cell", cells);
1520 0 : Vector<String> vcells(cells);
1521 0 : if(cells.nelements()==1) // [ cellx ]
1522 : {
1523 0 : Quantity qcell;
1524 0 : err+= stringToQuantity( vcells[0], qcell ); cellsize.set( qcell );
1525 : }
1526 0 : else if( cells.nelements()==2 ) // [ cellx, celly ]
1527 : {
1528 0 : err+= stringToQuantity( vcells[0], cellsize[0] );
1529 0 : 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 0 : err += readVal( inrec, String("stokes"), stokes);
1546 0 : if(stokes.matches("pseudoI"))
1547 : {
1548 0 : stokes="I";
1549 0 : pseudoi=true;
1550 : }
1551 0 : else {pseudoi=false;}
1552 :
1553 : /// PseudoI
1554 :
1555 : ////nchan
1556 0 : err += readVal( inrec, String("nchan"), nchan);
1557 :
1558 : /// phaseCenter (as a string) . // Add INT support later.
1559 : //err += readVal( inrec, String("phasecenter"), phaseCenter );
1560 0 : if( inrec.isDefined("phasecenter") )
1561 : {
1562 0 : String pcent("");
1563 0 : if( inrec.dataType("phasecenter") == TpString )
1564 : {
1565 0 : inrec.get("phasecenter",pcent);
1566 0 : if( pcent.length() > 0 ) // if it's zero length, it means 'figure out from first field in MS'.
1567 : {
1568 0 : err += readVal( inrec, String("phasecenter"), phaseCenter );
1569 0 : 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 0 : err += readVal( inrec, String("phasecenterfieldid"), phaseCenterFieldId);
1573 : }
1574 : //else { phaseCenterFieldId=0; } // Take the first field of the MS.
1575 0 : else { phaseCenterFieldId=-2; } // deal with this later in buildCoordinateSystem to assign the first selected field
1576 : }
1577 0 : if (inrec.dataType("phasecenter")==TpInt
1578 0 : || inrec.dataType("phasecenter")==TpFloat
1579 0 : || inrec.dataType("phasecenter")==TpDouble )
1580 : {
1581 : // This will override the previous setting to 0 if the phaseCenter string is zero length.
1582 0 : err += readVal( inrec, String("phasecenter"), phaseCenterFieldId );
1583 : }
1584 :
1585 0 : if( ( inrec.dataType("phasecenter") != TpString && inrec.dataType("phasecenter")!=TpInt
1586 0 : && 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 0 : if( inrec.isDefined("projection") )
1596 : {
1597 0 : if( inrec.dataType("projection") == TpString )
1598 : {
1599 0 : String pstr;
1600 0 : inrec.get("projection",pstr);
1601 :
1602 : try
1603 : {
1604 0 : if( pstr.matches("NCP") )
1605 : {
1606 0 : pstr ="SIN";
1607 0 : useNCP=true;
1608 : }
1609 0 : 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 0 : err += readVal( inrec, String("specmode"), mode);
1621 : // Alias for 'mfs' is 'cont'
1622 0 : if(mode=="cont") mode="mfs";
1623 :
1624 0 : err += readVal( inrec, String("outframe"), frame);
1625 0 : qmframe="";
1626 : // mveltype is only set when start/step is given in mdoppler
1627 0 : mveltype="";
1628 : //start
1629 0 : String startType("");
1630 0 : String widthType("");
1631 0 : if( inrec.isDefined("start") )
1632 : {
1633 0 : if( inrec.dataType("start") == TpInt )
1634 : {
1635 0 : err += readVal( inrec, String("start"), chanStart);
1636 0 : start = String::toString(chanStart);
1637 0 : startType="chan";
1638 : }
1639 0 : else if( inrec.dataType("start") == TpString )
1640 : {
1641 0 : err += readVal( inrec, String("start"), start);
1642 0 : if( start.contains("Hz") )
1643 : {
1644 0 : stringToQuantity(start,freqStart);
1645 0 : startType="freq";
1646 : }
1647 0 : else if( start.contains("m/s") )
1648 : {
1649 0 : stringToQuantity(start,velStart);
1650 0 : startType="vel";
1651 : }
1652 : }
1653 0 : 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 0 : startRecord = inrec.subRecord("start");
1659 0 : if(startRecord.isDefined("m0") )
1660 : {
1661 : //must be a measure
1662 0 : String mtype;
1663 0 : startRecord.get("type", mtype);
1664 0 : if( mtype=="frequency")
1665 : {
1666 : //mfrequency
1667 0 : startRecord.get(String("refer"), qmframe);
1668 0 : 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 0 : start = recordQMToString(startRecord);
1674 0 : stringToQuantity(start,freqStart);
1675 0 : startType="freq";
1676 : }
1677 0 : else if( mtype=="radialvelocity")
1678 : {
1679 : //mradialvelocity
1680 0 : startRecord.get(String("refer"), qmframe);
1681 0 : if ( frame!="" && frame!=qmframe)
1682 : {
1683 : // should emit warning to the logger
1684 0 : cerr<<"The frame in start:"<<qmframe<<" Override frame="<<frame<<endl;
1685 : }
1686 0 : start = recordQMToString(startRecord);
1687 0 : stringToQuantity(start,velStart);
1688 0 : startType="vel";
1689 : }
1690 0 : else if( mtype=="doppler")
1691 : {
1692 : //use veltype in mdoppler
1693 : //start = MDopToVelString(startRecord);
1694 0 : start = recordQMToString(startRecord);
1695 0 : stringToQuantity(start,velStart);
1696 0 : startRecord.get(String("refer"), mveltype);
1697 0 : mveltype.downcase();
1698 0 : startType="vel";
1699 : }
1700 : }
1701 : else
1702 : {
1703 0 : start = recordQMToString(startRecord);
1704 0 : if ( start.contains("Hz") )
1705 : {
1706 0 : stringToQuantity(start,freqStart);
1707 0 : startType="freq";
1708 : }
1709 0 : else if ( start.contains("m/s") )
1710 : {
1711 0 : stringToQuantity(start,velStart);
1712 0 : 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 0 : if( inrec.isDefined("width") )
1722 : {
1723 0 : if( inrec.dataType("width") == TpInt )
1724 : {
1725 0 : err += readVal( inrec, String("width"), chanStep);
1726 0 : step = String::toString(chanStep);
1727 0 : widthType="chan";
1728 : }
1729 0 : else if( inrec.dataType("width") == TpString )
1730 : {
1731 0 : err += readVal( inrec, String("width"), step);
1732 0 : if( step.contains("Hz") )
1733 : {
1734 0 : stringToQuantity(step,freqStep);
1735 0 : widthType="freq";
1736 : }
1737 0 : else if( step.contains("m/s") )
1738 : {
1739 0 : stringToQuantity(step,velStep);
1740 0 : widthType="vel";
1741 : }
1742 : }
1743 0 : 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 0 : stepRecord = inrec.subRecord("width");
1749 0 : if(stepRecord.isDefined("m0") )
1750 : {
1751 : //must be a measure
1752 0 : String mtype;
1753 0 : stepRecord.get("type", mtype);
1754 0 : 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 0 : else if( mtype=="radialvelocity")
1768 : {
1769 : //mradialvelocity
1770 0 : stepRecord.get(String("refer"), qmframe);
1771 0 : 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 0 : step = recordQMToString(stepRecord);
1777 0 : stringToQuantity(step,velStep);
1778 0 : 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 0 : step = recordQMToString(stepRecord);
1793 0 : if ( step.contains("Hz") )
1794 : {
1795 0 : stringToQuantity(step,freqStep);
1796 0 : widthType="freq";
1797 : }
1798 0 : else if ( step.contains("m/s") )
1799 : {
1800 0 : stringToQuantity(step,velStep);
1801 0 : 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 0 : 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 0 : if( inrec.isDefined("reffreq") )
1815 : {
1816 0 : if( inrec.dataType("reffreq")==TpString )
1817 : {
1818 0 : 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 0 : err += readVal( inrec, String("veltype"), veltype);
1846 0 : veltype = mveltype!=""? mveltype:veltype;
1847 : // sysvel (String, Quantity)
1848 0 : if( inrec.isDefined("sysvel") )
1849 : {
1850 0 : if( inrec.dataType("sysvel")==TpString )
1851 : {
1852 0 : 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 0 : err += readVal( inrec, String("sysvelframe"), sysvelframe);
1865 :
1866 : // rest frequencies (record or vector of Strings)
1867 0 : if( inrec.isDefined("restfreq") )
1868 : {
1869 0 : Vector<String> rfreqs(0);
1870 0 : Record restfreqSubRecord;
1871 0 : 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 0 : else if( inrec.dataType("restfreq")==TpArrayString )
1886 : {
1887 : //Vector<String> rfreqs(0);
1888 0 : err += readVal( inrec, String("restfreq"), rfreqs );
1889 : // case no restfreq is given: set to
1890 : }
1891 0 : else if( inrec.dataType("restfreq")==TpString )
1892 : {
1893 0 : rfreqs.resize(1);
1894 0 : err += readVal( inrec, String("restfreq"), rfreqs[0] );
1895 : // case no restfreq is given: set to
1896 : }
1897 0 : restFreq.resize( rfreqs.nelements() );
1898 0 : for( uInt fr=0; fr<rfreqs.nelements(); fr++)
1899 : {
1900 0 : 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 0 : if( inrec.isDefined("csys") )
1907 : {
1908 : // cout<<"HAS CSYS KEY - got from input record"<<endl;
1909 0 : if( inrec.dataType("csys")==TpRecord )
1910 : {
1911 : //csysRecord = inrec.subRecord("csys");
1912 0 : csysRecord.defineRecord("coordsys",inrec.subRecord("csys"));
1913 : }
1914 0 : if( inrec.isDefined("imshape") )
1915 : {
1916 0 : if ( inrec.dataType("imshape") == TpArrayInt )
1917 : {
1918 0 : 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 0 : String freqframestr = (frame!="" && qmframe!="")? qmframe:frame;
1929 0 : if( frame!="" && ! MFrequency::getType(freqFrame, freqframestr) )
1930 0 : { err += "Invalid Frequency Frame " + freqframestr ; }
1931 0 : err += readVal( inrec, String("restart"), overwrite );
1932 :
1933 0 : err += readVal(inrec, String("freqframevalid"), freqFrameValid);
1934 : // startmodel parsing copied in SynthesisParamDeconv. Clean this up !!!
1935 0 : if( inrec.isDefined("startmodel") )
1936 : {
1937 0 : if( inrec.dataType("startmodel")==TpString )
1938 : {
1939 0 : String onemodel;
1940 0 : err += readVal( inrec, String("startmodel"), onemodel );
1941 0 : if( onemodel.length()>0 )
1942 : {
1943 0 : startModel.resize(1);
1944 0 : startModel[0] = onemodel;
1945 : }
1946 0 : else {startModel.resize();}
1947 : }
1948 0 : else if( inrec.dataType("startmodel")==TpArrayString ||
1949 0 : inrec.dataType("startmodel")==TpArrayBool)
1950 : {
1951 0 : 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 0 : err += readVal( inrec, String("nterms"), nTaylorTerms );
1959 0 : err += readVal( inrec, String("deconvolver"), deconvolver );
1960 :
1961 : // Force nchan=1 for anything other than cube modes...
1962 0 : if(mode=="mfs") nchan=1;
1963 : //read obslocation
1964 0 : if(inrec.isDefined("obslocation_rec")){
1965 0 : String errorobs;
1966 0 : const Record obsrec=inrec.asRecord("obslocation_rec");
1967 0 : MeasureHolder mh;
1968 0 : if(!mh.fromRecord(errorobs, obsrec)){
1969 0 : err+=errorobs;
1970 : }
1971 0 : obslocation=mh.asMPosition();
1972 :
1973 : }
1974 :
1975 :
1976 :
1977 0 : err += verify();
1978 :
1979 : }
1980 0 : catch(AipsError &x)
1981 : {
1982 0 : err = err + x.getMesg() + "\n";
1983 : }
1984 :
1985 0 : if( err.length()>0 ) throw(AipsError("Invalid Image Parameter set : " + err));
1986 :
1987 0 : }
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 0 : String SynthesisParamsImage::verify() const
2021 : {
2022 0 : String err;
2023 :
2024 0 : if( imageName=="" ) {err += "Please supply an image name\n";}
2025 :
2026 0 : if( imsize.nelements() != 2 ){ err += "imsize must be a vector of 2 Ints\n"; }
2027 0 : if( cellsize.nelements() != 2 ) { err += "cellsize must be a vector of 2 Quantities\n"; }
2028 0 : if( cellsize[0].getValue() == 0.0 || cellsize[1].getValue() == 0.0 ) {
2029 0 : 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 0 : if( (mode=="mfs") && nchan>1 )
2037 0 : { err += "specmode=mfs cannot have nchan="+String::toString(nchan)+" (must be 1)\n";}
2038 :
2039 0 : if( ! stokes.matches("I") && ! stokes.matches("Q") &&
2040 0 : ! stokes.matches("U") && ! stokes.matches("V") &&
2041 0 : ! stokes.matches("RR") && ! stokes.matches("LL") &&
2042 0 : ! stokes.matches("XX") && ! stokes.matches("YY") &&
2043 0 : ! stokes.matches("IV") && ! stokes.matches("IQ") &&
2044 0 : ! stokes.matches("RRLL") && ! stokes.matches("XXYY") &&
2045 0 : ! stokes.matches("QU") && ! stokes.matches("UV") &&
2046 0 : ! stokes.matches("IQU") && ! stokes.matches("IUV") &&
2047 0 : ! 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 0 : if( startModel.nelements()>0 )
2054 : {
2055 0 : if( deconvolver!="mtmfs" ) {
2056 :
2057 0 : if( startModel.nelements()!=1 ){err += String("Only one startmodel image is allowed.\n");}
2058 : else
2059 : {
2060 0 : File fp( imageName+String(".model") );
2061 0 : if( fp.exists() ) err += "Model " + imageName+".model exists, but a starting model of " + startModel[0] + " is also being requested. Please either reset startmodel='' to use what already exists, or delete " + imageName + ".model so that it uses the new model specified in startmodel.";
2062 : }
2063 : }
2064 : else {// mtmfs
2065 0 : File fp( imageName+String(".model.tt0") );
2066 0 : 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 0 : for(uInt ss=0;ss<startModel.nelements();ss++)
2075 : {
2076 0 : File fp( startModel[ss] );
2077 0 : if( ! fp.exists() ) {err += "Startmodel " + startModel[ss] + " cannot be found on disk.";}
2078 : }
2079 :
2080 : }
2081 :
2082 :
2083 : /// Check imsize for efficiency.
2084 0 : Int imxnew = SynthesisUtilMethods::getOptimumSize( imsize[0] );
2085 0 : Int imynew = SynthesisUtilMethods::getOptimumSize( imsize[1] );
2086 :
2087 0 : if( imxnew != imsize[0] || imynew != imsize[1] )
2088 : {
2089 0 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2090 0 : if( imxnew != imsize[0] ) {os << LogIO::WARN << "imsize with "+String::toString(imsize[0])+" pixels is not an efficient imagesize. Try "+String::toString(imxnew)+" instead." << LogIO::POST;}
2091 0 : if( imsize[0] != imsize[1] && imynew != imsize[1] ) {os << LogIO::WARN << "imsize with "+String::toString(imsize[1])+" pixels is not an efficient imagesize. Try "+String::toString(imynew)+" instead." << LogIO::POST;}
2092 : //err += "blah";
2093 : }
2094 :
2095 0 : 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 0 : void SynthesisParamsImage::setDefaults()
2109 : {
2110 : // Image definition parameters
2111 0 : imageName = String("");
2112 0 : imsize.resize(2); imsize.set(100);
2113 0 : cellsize.resize(2); cellsize.set( Quantity(1.0,"arcsec") );
2114 0 : stokes="I";
2115 0 : phaseCenter=MDirection();
2116 0 : phaseCenterFieldId=-1;
2117 0 : projection=Projection::SIN;
2118 0 : useNCP=false;
2119 0 : startModel=Vector<String>(0);
2120 0 : freqFrameValid=True;
2121 0 : overwrite=false;
2122 : // PseudoI
2123 0 : pseudoi=false;
2124 :
2125 : // Spectral coordinates
2126 0 : nchan=1;
2127 0 : mode="mfs";
2128 0 : start="";
2129 0 : step="";
2130 0 : chanStart=0;
2131 0 : 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 0 : freqStart=Quantity(0,"");
2137 0 : freqStep=Quantity(0,"");
2138 0 : velStart=Quantity(0,"");
2139 0 : velStep=Quantity(0,"");
2140 0 : veltype=String("radio");
2141 0 : restFreq.resize(0);
2142 0 : refFreq = Quantity(0,"Hz");
2143 0 : frame = "";
2144 0 : freqFrame=MFrequency::LSRK;
2145 0 : sysvel="";
2146 0 : sysvelframe="";
2147 0 : sysvelvalue=Quantity(0.0,"m/s");
2148 0 : nTaylorTerms=1;
2149 0 : deconvolver="hogbom";
2150 : ///csysRecord=Record();
2151 : //
2152 :
2153 :
2154 0 : }
2155 :
2156 0 : Record SynthesisParamsImage::toRecord() const
2157 : {
2158 0 : Record impar;
2159 0 : impar.define("imagename", imageName);
2160 0 : impar.define("imsize", imsize);
2161 0 : Vector<String> cells(2);
2162 0 : cells[0] = QuantityToString( cellsize[0] );
2163 0 : cells[1] = QuantityToString( cellsize[1] );
2164 0 : impar.define("cell", cells );
2165 0 : if(pseudoi==true){impar.define("stokes","pseudoI");}
2166 0 : else{impar.define("stokes", stokes);}
2167 0 : impar.define("nchan", nchan);
2168 0 : impar.define("nterms", nTaylorTerms);
2169 0 : impar.define("deconvolver",deconvolver);
2170 0 : impar.define("phasecenter", MDirectionToString( phaseCenter ) );
2171 0 : impar.define("phasecenterfieldid",phaseCenterFieldId);
2172 0 : impar.define("projection", (useNCP? "NCP" : projection.name()) );
2173 :
2174 0 : impar.define("specmode", mode );
2175 : // start and step can be one of these types
2176 0 : if( start!="" )
2177 : {
2178 0 : if( !start.contains("Hz") && !start.contains("m/s") &&
2179 0 : String::toInt(start) == chanStart )
2180 : {
2181 0 : impar.define("start",chanStart);
2182 : }
2183 0 : else if( startRecord.nfields() > 0 )
2184 : {
2185 0 : impar.defineRecord("start", startRecord );
2186 : }
2187 : else
2188 : {
2189 0 : impar.define("start",start);
2190 : }
2191 : }
2192 : else {
2193 0 : impar.define("start", start );
2194 : }
2195 0 : if( step!="" )
2196 : {
2197 0 : if( !step.contains("Hz") && !step.contains("m/s") &&
2198 0 : String::toInt(step) == chanStep )
2199 : {
2200 0 : impar.define("width", chanStep);
2201 : }
2202 0 : else if( stepRecord.nfields() > 0 )
2203 : {
2204 0 : impar.defineRecord("width",stepRecord);
2205 : }
2206 : else
2207 : {
2208 0 : impar.define("width",step);
2209 : }
2210 : }
2211 : else
2212 : {
2213 0 : 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 0 : impar.define("veltype", veltype);
2222 0 : if (restfreqRecord.nfields() != 0 )
2223 : {
2224 0 : impar.defineRecord("restfreq", restfreqRecord);
2225 : }
2226 : else
2227 : {
2228 0 : Vector<String> rfs( restFreq.nelements() );
2229 0 : for(uInt rf=0; rf<restFreq.nelements(); rf++){rfs[rf] = QuantityToString(restFreq[rf]);}
2230 0 : impar.define("restfreq", rfs);
2231 : }
2232 : //impar.define("reffreq", QuantityToString(refFreq));
2233 : //reffreq
2234 0 : if( reffreqRecord.nfields() != 0 )
2235 0 : { impar.defineRecord("reffreq",reffreqRecord); }
2236 : else
2237 0 : { impar.define("reffreq",reffreq); }
2238 : //impar.define("reffreq", reffreq );
2239 : //impar.define("outframe", MFrequency::showType(freqFrame) );
2240 0 : impar.define("outframe", frame );
2241 : //sysvel
2242 0 : if( sysvelRecord.nfields() != 0 )
2243 0 : { impar.defineRecord("sysvel",sysvelRecord); }
2244 : else
2245 0 : { impar.define("sysvel", sysvel );}
2246 0 : impar.define("sysvelframe", sysvelframe );
2247 :
2248 0 : impar.define("restart",overwrite );
2249 0 : impar.define("freqframevalid", freqFrameValid);
2250 0 : impar.define("startmodel", startModel );
2251 :
2252 0 : if( csysRecord.isDefined("coordsys") )
2253 : {
2254 : // cout <<" HAS CSYS INFO.... writing to output record"<<endl;
2255 0 : impar.defineRecord("csys", csysRecord.subRecord("coordsys"));
2256 0 : impar.define("imshape", imshape);
2257 : }
2258 : // else cout << " NO CSYS INFO to write to output record " << endl;
2259 : ///Now save obslocation
2260 0 : Record tmprec;
2261 0 : String err;
2262 0 : MeasureHolder mh(obslocation);
2263 0 : if(mh.toRecord(err, tmprec)){
2264 0 : impar.defineRecord("obslocation_rec", tmprec);
2265 : }
2266 : else{
2267 0 : throw(AipsError("failed to save obslocation to record"));
2268 : }
2269 0 : 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 0 : 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 0 : vi::VisBuffer2* vb=vi2.getVisBuffer();
2290 0 : vi2.originChunks();
2291 0 : vi2.origin();
2292 : /// This version uses the new vi2/vb2
2293 : // get the first ms for multiple MSes
2294 : //MeasurementSet msobj=vi2.ms();
2295 0 : Int fld=vb->fieldId()(0);
2296 :
2297 : //handling first ms only
2298 0 : Double gfreqmax=-1.0;
2299 0 : Double gdatafend=-1.0;
2300 0 : Double gdatafstart=1e14;
2301 0 : Double gfreqmin=1e14;
2302 0 : Vector<Int> spwids0;
2303 0 : Int j=0;
2304 0 : Int minfmsid=0;
2305 : //for cube mode ,for a list of MSs, check ms to send to buildCoordSysCore contains start freq/vel
2306 0 : Double imStartFreq=getCubeImageStartFreq();
2307 0 : std::vector<Int> sourceMsWithStartFreq;
2308 :
2309 :
2310 0 : for (auto forMS0=chansel.begin(); forMS0 !=chansel.end(); ++forMS0, ++j){
2311 : //auto forMS0=chansel.find(0);
2312 0 : map<Int, Vector<Int> > spwsels=forMS0->second;
2313 0 : Int nspws=spwsels.size();
2314 0 : Vector<Int> spwids(nspws);
2315 0 : Vector<Int> nChannels(nspws);
2316 0 : Vector<Int> firstChannels(nspws);
2317 : //Vector<Int> channelIncrement(nspws);
2318 :
2319 0 : Int k=0;
2320 0 : for (auto it=spwsels.begin(); it != spwsels.end(); ++it, ++k){
2321 0 : spwids[k]=it->first;
2322 0 : nChannels[k]=(it->second)[0];
2323 0 : firstChannels[k]=(it->second)[1];
2324 : }
2325 0 : if(j==0) {
2326 0 : spwids0.resize();
2327 0 : 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 0 : Double freqmin=0, freqmax=0;
2339 0 : freqFrameValid=(freqFrame != MFrequency::REST || mode=="cubesource");
2340 :
2341 : //MFrequency::Types dataFrame=(MFrequency::Types)vi2.subtableColumns().spectralWindow().measFreqRef()(spwids[0]);
2342 0 : 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 0 : Bool status=MSUtil::getFreqRangeInSpw( datafstart, datafend, spwids, firstChannels, nChannels,*mss[j], dataFrame, True);
2348 : //cerr << "after " << datafstart << " " << datafend << endl;
2349 0 : if((datafstart > datafend) || !status)
2350 0 : throw(AipsError("spw selection failed"));
2351 : //cerr << "datafstart " << datafstart << " end " << datafend << endl;
2352 :
2353 0 : if (mode=="cubedata") {
2354 0 : freqmin = datafstart;
2355 0 : freqmax = datafend;
2356 : }
2357 0 : else if(mode == "cubesource"){
2358 0 : if(!trackSource){
2359 0 : throw(AipsError("Cannot be in cubesource without tracking a moving source"));
2360 : }
2361 0 : String ephemtab(movingSource);
2362 0 : if(movingSource=="TRACKFIELD"){
2363 0 : Int fieldID=MSColumns(*mss[j]).fieldId()(0);
2364 0 : ephemtab=Path(MSColumns(*mss[j]).field().ephemPath(fieldID)).absoluteName();
2365 : }
2366 0 : MEpoch refep=MSColumns(*mss[j]).timeMeas()(0);
2367 0 : Quantity refsysvel;
2368 0 : MSUtil::getFreqRangeAndRefFreqShift(freqmin,freqmax,refsysvel, refep, spwids,firstChannels, nChannels, *mss[j], ephemtab, trackDir, true);
2369 0 : if(j==0)
2370 0 : 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 0 : MSUtil::getFreqRangeInSpw( freqmin, freqmax, spwids, firstChannels,
2382 0 : nChannels,*mss[j], freqFrameValid? freqFrame:MFrequency::REST , True);
2383 : //cerr << "after " << freqmin << " " << freqmax << endl;
2384 : }
2385 :
2386 :
2387 :
2388 :
2389 0 : if(freqmin < gfreqmin) gfreqmin=freqmin;
2390 0 : if(freqmax > gfreqmax) gfreqmax=freqmax;
2391 0 : if(datafstart < gdatafstart) gdatafstart=datafstart;
2392 0 : 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 0 : if(imStartFreq > 0.0 && imStartFreq >= freqmin && imStartFreq <= freqmax){
2399 0 : if(mode != "cubesource"){
2400 0 : minfmsid=j;
2401 0 : spwids0.resize();
2402 0 : spwids0=spwids;
2403 0 : vi2.originChunks();
2404 0 : vi2.origin();
2405 0 : while(vb->msId() != j && vi2.moreChunks() ){
2406 0 : vi2.nextChunk();
2407 0 : vi2.origin();
2408 : }
2409 0 : fld=vb->fieldId()(0);
2410 :
2411 : }
2412 : else{
2413 0 : sourceMsWithStartFreq.push_back(j);
2414 : }
2415 : }
2416 :
2417 : }
2418 0 : 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 0 : MeasurementSet msobj = *mss[minfmsid];
2425 : // return buildCoordinateSystemCore( msobj, spwids0, fld, gfreqmin, gfreqmax, gdatafstart, gdatafend );
2426 0 : 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 0 : CoordinateSystem SynthesisParamsImage::buildCoordinateSystemCore(
2460 : MeasurementSet& msobj,
2461 : Vector<Int> spwids, Int fld,
2462 : Double freqmin, Double freqmax,
2463 : Double datafstart, Double datafend )
2464 : {
2465 0 : LogIO os( LogOrigin("SynthesisParamsImage","buildCoordinateSystem",WHERE) );
2466 :
2467 0 : CoordinateSystem csys;
2468 0 : if( csysRecord.nfields()!=0 )
2469 : {
2470 : //use cysRecord
2471 0 : Record subRec1;
2472 : // cout<<"USE THE EXISTING CSYS +++++++++++++++++"<<endl;
2473 0 : CoordinateSystem *csysptr = CoordinateSystem::restore(csysRecord,"coordsys");
2474 : //csys = *csysptr;
2475 : //CoordinateSystem csys(*csysptr);
2476 0 : csys = *csysptr;
2477 :
2478 : }
2479 : else {
2480 0 : MSColumns msc(msobj);
2481 0 : String telescop = msc.observation().telescopeName()(0);
2482 0 : MEpoch obsEpoch = msc.timeMeas()(0);
2483 0 : MPosition obsPosition;
2484 0 : 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 0 : MDirection phaseCenterToUse = phaseCenter;
2496 :
2497 0 : if( phaseCenterFieldId != -1 )
2498 : {
2499 0 : MSFieldColumns msfield(msobj.field());
2500 0 : if(phaseCenterFieldId == -2) // the case for phasecenter=''
2501 : {
2502 0 : if(trackSource){
2503 0 : phaseCenterToUse=getMovingSourceDir(msobj, obsEpoch, obsPosition, MDirection::ICRS);
2504 : }
2505 : else{
2506 0 : phaseCenterToUse=msfield.phaseDirMeas( fld );
2507 : }
2508 : }
2509 : else
2510 : {
2511 0 : phaseCenterToUse=msfield.phaseDirMeas( phaseCenterFieldId );
2512 : }
2513 : }
2514 : // Setup Phase center if it is specified only by field id.
2515 :
2516 : /////////////////// Direction Coordinates
2517 0 : MVDirection mvPhaseCenter(phaseCenterToUse.getAngle());
2518 : // Normalize correctly
2519 0 : MVAngle ra=mvPhaseCenter.get()(0);
2520 0 : ra(0.0);
2521 :
2522 0 : MVAngle dec=mvPhaseCenter.get()(1);
2523 0 : Vector<Double> refCoord(2);
2524 0 : refCoord(0)=ra.get().getValue();
2525 0 : refCoord(1)=dec;
2526 0 : Vector<Double> refPixel(2);
2527 0 : refPixel(0) = Double(imsize[0]/2);
2528 0 : refPixel(1) = Double(imsize[1]/2);
2529 :
2530 0 : Vector<Double> deltas(2);
2531 0 : deltas(0)=-1* cellsize[0].get("rad").getValue();
2532 0 : deltas(1)=cellsize[1].get("rad").getValue();
2533 0 : Matrix<Double> xform(2,2);
2534 0 : xform=0.0;xform.diagonal()=1.0;
2535 :
2536 0 : Vector<Double> projparams(2);
2537 0 : projparams = 0.0;
2538 0 : if( useNCP==true ) { projparams[0]=0.0, projparams[1]=1/tan(refCoord(1)); }
2539 0 : Projection projTo( projection.type(), projparams );
2540 :
2541 : DirectionCoordinate
2542 0 : myRaDec(MDirection::Types(phaseCenterToUse.getRefPtr()->getType()),
2543 : // projection,
2544 : projTo,
2545 0 : refCoord(0), refCoord(1),
2546 0 : deltas(0), deltas(1),
2547 : xform,
2548 0 : 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 0 : obslocation=obsPosition;
2556 0 : ObsInfo myobsinfo;
2557 0 : myobsinfo.setTelescope(telescop);
2558 0 : myobsinfo.setPointingCenter(mvPhaseCenter);
2559 0 : myobsinfo.setObsDate(obsEpoch);
2560 0 : 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 0 : if(spwids.nelements()==0)
2576 : {
2577 0 : Int nspw=msc.spectralWindow().nrow();
2578 0 : spwids.resize(nspw);
2579 0 : indgen(spwids);
2580 : }
2581 0 : MFrequency::Types dataFrame=(MFrequency::Types)msc.spectralWindow().measFreqRef()(spwids[0]);
2582 0 : Vector<Double> dataChanFreq, dataChanWidth;
2583 0 : std::vector<std::vector<Int> > averageWhichChan;
2584 0 : std::vector<std::vector<Int> > averageWhichSPW;
2585 0 : std::vector<std::vector<Double> > averageChanFrac;
2586 :
2587 0 : if(spwids.nelements()==1)
2588 : {
2589 0 : dataChanFreq=msc.spectralWindow().chanFreq()(spwids[0]);
2590 0 : dataChanWidth=msc.spectralWindow().chanWidth()(spwids[0]);
2591 : }
2592 : else
2593 : {
2594 0 : 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 0 : Double minDataFreq = min(dataChanFreq);
2601 0 : if(start=="" && minDataFreq < datafstart ) {
2602 : // limit data chan freq vector for default start case with channel selection
2603 : Int chanStart, chanEnd;
2604 0 : Int lochan = 0;
2605 0 : Int nDataChan = dataChanFreq.nelements();
2606 0 : Int hichan = nDataChan-1;
2607 : Double diff_fmin, diff_fmax;
2608 0 : Bool ascending = dataChanFreq[nDataChan-1] - dataChanFreq[0] > 0;
2609 0 : for(Int ichan = 0; ichan < nDataChan; ichan++)
2610 : {
2611 0 : diff_fmin = dataChanFreq[ichan] - datafstart;
2612 0 : diff_fmax = datafend - dataChanFreq[ichan];
2613 : // freqmin and freqmax should corresponds to the channel edges
2614 0 : if(ascending)
2615 : {
2616 :
2617 0 : if( diff_fmin > 0 && diff_fmin <= dataChanWidth[ichan]/2. )
2618 : {
2619 0 : lochan = ichan;
2620 : }
2621 0 : else if(diff_fmax > 0 && diff_fmax <= dataChanWidth[ichan]/2. )
2622 : {
2623 0 : 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 0 : chanStart = lochan;
2639 0 : chanEnd = hichan;
2640 0 : if (lochan > hichan)
2641 : {
2642 0 : chanStart=hichan;
2643 0 : chanEnd=lochan;
2644 : }
2645 0 : Vector<Double> tempChanFreq = dataChanFreq(Slice(chanStart,chanEnd-chanStart+1,1));
2646 0 : Vector<Double> tempChanWidth = dataChanWidth(Slice(chanStart,chanEnd-chanStart+1,1));
2647 0 : dataChanFreq.resize(tempChanFreq.nelements());
2648 0 : dataChanWidth.resize(tempChanWidth.nelements());
2649 0 : dataChanFreq = tempChanFreq;
2650 0 : dataChanWidth = tempChanWidth;
2651 : }
2652 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
2653 0 : String cubemode;
2654 0 : if ( qrestfreq.getValue("Hz")==0 )
2655 : {
2656 0 : MSDopplerUtil msdoppler(msobj);
2657 0 : Vector<Double> restfreqvec;
2658 0 : msdoppler.dopplerInfo(restfreqvec, spwids(0), fld);
2659 0 : qrestfreq = restfreqvec.nelements() >0 ? Quantity(restfreqvec(0),"Hz"): Quantity(0.0, "Hz");
2660 0 : if ( qrestfreq.getValue("Hz")==0 and mode!="mfs" )
2661 : {
2662 0 : cubemode = findSpecMode(mode);
2663 0 : 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 0 : Double provisional_restfreq = (datafend+datafstart)/2.0;
2668 0 : 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 0 : << 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 0 : Vector<Double> chanFreq;
2680 0 : Vector<Double> chanFreqStep;
2681 0 : String specmode;
2682 :
2683 0 : if(mode=="cubesource"){
2684 0 : MDoppler mdop(sysvelvalue, MDoppler::RELATIVISTIC);
2685 0 : dataChanFreq=mdop.shiftFrequency(dataChanFreq);
2686 0 : dataChanWidth=mdop.shiftFrequency(dataChanWidth);
2687 0 : 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 0 : 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 0 : Bool nonLinearFreq(false);
2699 0 : String veltype_p=veltype;
2700 0 : veltype_p.upcase();
2701 0 : if(veltype_p.contains("OPTICAL") || veltype_p.matches("Z") || veltype_p.contains("BETA") ||
2702 0 : veltype_p.contains("RELATI") || veltype_p.contains("GAMMA"))
2703 : {
2704 0 : nonLinearFreq= true;
2705 : }
2706 :
2707 0 : SpectralCoordinate mySpectral;
2708 : Double stepf;
2709 0 : if(!nonLinearFreq)
2710 : //TODO: velocity mode default start case (use last channels?)
2711 : {
2712 0 : Double startf=chanFreq[0];
2713 : //Double stepf=chanFreqStep[0];
2714 0 : if(chanFreq.nelements()==1)
2715 : {
2716 0 : stepf=chanFreqStep[0];
2717 : }
2718 : else
2719 : {
2720 0 : stepf=chanFreq[1]-chanFreq[0];
2721 : }
2722 0 : Double restf=qrestfreq.getValue("Hz");
2723 : //stepf=9e8;
2724 0 : 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 0 : if(mode=="cubedata")
2728 : {
2729 : // mySpectral = SpectralCoordinate(freqFrameValid ? MFrequency::Undefined : MFrequency::REST,
2730 0 : mySpectral = SpectralCoordinate(freqFrame == MFrequency::REST?
2731 : MFrequency::REST : MFrequency::Undefined,
2732 0 : startf, stepf, refPix, restf);
2733 : }
2734 0 : 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 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
2740 0 : startf, stepf, refPix, restf);
2741 : }
2742 : else
2743 : {
2744 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
2745 0 : startf, stepf, refPix, restf);
2746 : }
2747 : }
2748 : else
2749 : { // nonlinear freq coords - use tabular setting
2750 : // once NOFRAME is implemented do this
2751 0 : 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 0 : else if (mode=="cubesource")
2759 : {
2760 0 : mySpectral = SpectralCoordinate(MFrequency::REST,
2761 0 : chanFreq, (Double)qrestfreq.getValue("Hz"));
2762 : }
2763 : else
2764 : {
2765 0 : mySpectral = SpectralCoordinate(freqFrameValid ? freqFrame : MFrequency::REST,
2766 0 : 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 0 : uInt nrestfreq = restFreq.nelements();
2775 0 : 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 0 : Vector<Int> whichStokes = decideNPolPlanes(stokes);
2793 0 : if(whichStokes.nelements()==0)
2794 0 : throw(AipsError("Stokes selection of " + stokes + " is invalid"));
2795 0 : StokesCoordinate myStokes(whichStokes);
2796 :
2797 : ////////////////// Build Full coordinate system.
2798 :
2799 : //CoordinateSystem csys;
2800 0 : csys.addCoordinate(myRaDec);
2801 0 : csys.addCoordinate(myStokes);
2802 0 : csys.addCoordinate(mySpectral);
2803 0 : csys.setObsInfo(myobsinfo);
2804 :
2805 : //store back csys to impars record
2806 : //cerr<<"save csys to csysRecord..."<<endl;
2807 0 : if(csysRecord.isDefined("coordsys"))
2808 0 : csysRecord.removeField("coordsys");
2809 0 : csys.save(csysRecord,"coordsys");
2810 : //cerr<<"BUILDCOORDSYS:: new csysRecord ="<<csysRecord<<endl;
2811 : // imshape
2812 0 : imshape.resize(4);
2813 0 : imshape[0] = imsize[0];
2814 0 : imshape[1] = imsize[0];
2815 0 : imshape[2] = whichStokes.nelements();
2816 0 : 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 0 : 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 0 : MDirection SynthesisParamsImage::getMovingSourceDir(const MeasurementSet& ms, const MEpoch& refEp, const MPosition& obsposition, const MDirection::Types outframe){
3106 0 : MDirection outdir;
3107 0 : String ephemtab(movingSource);
3108 0 : if(movingSource=="TRACKFIELD"){
3109 0 : Int fieldID=MSColumns(ms).fieldId()(0);
3110 0 : ephemtab=Path(MSColumns(ms).field().ephemPath(fieldID)).absoluteName();
3111 : }
3112 0 : casacore::MDirection::Types planetType=MDirection::castType(trackDir.getRef().getType());
3113 0 : 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 0 : MeasFrame mframe(refEp, obsposition);
3116 0 : MDirection::Ref outref1(MDirection::AZEL, mframe);
3117 0 : MDirection::Ref outref(outframe, mframe);
3118 0 : MDirection tmpazel;
3119 0 : if(planetType >=MDirection::MERCURY && planetType <MDirection::COMET){
3120 0 : tmpazel=MDirection::Convert(trackDir, outref1)();
3121 : }
3122 : else{
3123 0 : MeasComet mcomet(Path(ephemtab).absoluteName());
3124 0 : mframe.set(mcomet);
3125 0 : tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
3126 : }
3127 0 : outdir=MDirection::Convert(tmpazel, outref)();
3128 :
3129 0 : return outdir;
3130 : }
3131 :
3132 0 : 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 0 : String inStart, inStep;
3143 0 : specmode = findSpecMode(mode);
3144 0 : String freqframe;
3145 0 : Bool verbose("true"); // verbose logging messages from calcChanFreqs
3146 0 : LogIO os( LogOrigin("SynthesisParamsImage","getImFreq",WHERE) );
3147 :
3148 0 : refPix=0.0;
3149 0 : Bool descendingfreq(false);
3150 0 : Bool descendingoutfreq(false);
3151 :
3152 0 : if( mode.contains("cube") )
3153 : {
3154 0 : String restfreq=QuantityToString(qrestfreq);
3155 : // use frame from input start or width in MFreaquency or MRadialVelocity
3156 0 : freqframe = qmframe!=""? qmframe: MFrequency::showType(freqFrame);
3157 : // emit warning here if qmframe is used
3158 : //
3159 0 : inStart = start;
3160 0 : inStep = step;
3161 0 : if( specmode=="channel" )
3162 : {
3163 0 : inStart = String::toString(chanStart);
3164 0 : inStep = String::toString(chanStep);
3165 : // negative step -> descending channel indices
3166 0 : if (inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3167 : // input frame is the data frame
3168 : //freqframe = MFrequency::showType(dataFrame);
3169 : }
3170 0 : 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 0 : if(inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3181 : }
3182 0 : 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 0 : if (!inStep.contains(casacore::Regex("^-"))) descendingfreq=true;
3194 : }
3195 :
3196 0 : if (inStep=='0') inStep="";
3197 :
3198 0 : MRadialVelocity mSysVel;
3199 0 : Quantity qVel;
3200 : MRadialVelocity::Types mRef;
3201 0 : if(mode!="cubesource")
3202 : {
3203 :
3204 :
3205 0 : 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 0 : freqframe=MFrequency::showType(dataFrame);
3215 0 : 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 0 : { mSysVel=MRadialVelocity();}
3222 : }
3223 : // cubedata mode: input start, step are those of the input data frame
3224 0 : if ( mode=="cubedata" )
3225 : {
3226 0 : freqframe=MFrequency::showType(dataFrame);
3227 0 : 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 0 : os << LogIO::DEBUG1<<"mode="<<mode<<" specmode="<<specmode<<" inStart="<<inStart
3238 : <<" inStep="<<inStep<<" restfreq="<<restfreq<<" freqframe="<<freqframe
3239 0 : <<" dataFrame="<<dataFrame <<" veltype="<<veltype<<" nchan="<<nchan
3240 0 : << LogIO::POST;
3241 0 : ostringstream ostr;
3242 0 : ostr << " phaseCenter='" << phaseCenter;
3243 0 : os << String(ostr)<<"' ";
3244 :
3245 : Double dummy; // dummy variable - weightScale is not used here
3246 0 : 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 0 : veltype,
3263 : verbose,
3264 : mSysVel
3265 : );
3266 :
3267 0 : if( nchan==-1 )
3268 : {
3269 0 : nchan = chanFreq.nelements();
3270 0 : os << LogIO::DEBUG1 << "Setting nchan to number of selected channels : " << nchan << LogIO::POST;
3271 : }
3272 :
3273 0 : if (!rst) {
3274 : os << LogIO::SEVERE << "calcChanFreqs failed, check input start and width parameters"
3275 0 : << LogIO::EXCEPTION;
3276 0 : return false;
3277 : }
3278 : os << LogIO::DEBUG1
3279 0 : <<"chanFreq 0="<<chanFreq[0]<<" chanFreq last="<<chanFreq[chanFreq.nelements()-1]
3280 0 : << LogIO::POST;
3281 :
3282 0 : if (chanFreq[0]>chanFreq[chanFreq.nelements()-1]) {
3283 0 : descendingoutfreq = true;
3284 : }
3285 :
3286 : //if (descendingfreq && !descendingoutfreq) {
3287 0 : if ((specmode=="channel" && descendingfreq==1)
3288 0 : || (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 0 : std::vector<Double> stlchanfreq;
3294 0 : chanFreq.tovector(stlchanfreq);
3295 0 : std::reverse(stlchanfreq.begin(),stlchanfreq.end());
3296 0 : chanFreq=Vector<Double>(stlchanfreq);
3297 0 : chanFreqStep=-chanFreqStep;
3298 : }
3299 : }
3300 0 : else if ( mode=="mfs" ) {
3301 0 : chanFreq.resize(1);
3302 0 : chanFreqStep.resize(1);
3303 : //chanFreqStep[0] = freqmax - freqmin;
3304 0 : Double freqmean = (freqmin + freqmax)/2;
3305 0 : if (refFreq.getValue("Hz")==0) {
3306 0 : chanFreq[0] = freqmean;
3307 0 : refPix = 0.0;
3308 0 : chanFreqStep[0] = freqmax - freqmin;
3309 : }
3310 : else {
3311 0 : chanFreq[0] = refFreq.getValue("Hz");
3312 : // Set the new reffreq to be the refPix (CAS-9518)
3313 0 : refPix = 0.0; // (refFreq.getValue("Hz") - freqmean)/chanFreqStep[0];
3314 : // A larger bandwidth to compensate for the shifted reffreq (CAS-9518)
3315 0 : chanFreqStep[0] = freqmax - freqmin + 2*fabs(chanFreq[0] - freqmean);
3316 : }
3317 :
3318 0 : if( nchan==-1 ) nchan=1;
3319 0 : if( qrestfreq.getValue("Hz")==0.0 ) {
3320 0 : restFreq.resize(1);
3321 0 : 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 0 : return true;
3331 :
3332 : }//getImFreq
3333 : /////////////////////////
3334 0 : Double SynthesisParamsImage::getCubeImageStartFreq(){
3335 0 : Double inStartFreq=-1.0;
3336 0 : String checkspecmode("");
3337 0 : if(mode.contains("cube")) {
3338 0 : checkspecmode = findSpecMode(mode);
3339 : }
3340 0 : if(checkspecmode!="") {
3341 0 : MFrequency::Types mfreqframe = frame!="" ? MFrequency::typeFromString(frame):MFrequency::LSRK;
3342 0 : if(checkspecmode=="channel") {
3343 0 : inStartFreq=-1.0;
3344 : }
3345 : else {
3346 0 : if(checkspecmode=="frequency") {
3347 0 : inStartFreq = freqStart.get("Hz").getValue();
3348 : }
3349 0 : else if(checkspecmode=="velocity") {
3350 : MDoppler::Types DopType;
3351 0 : MDoppler::getType(DopType, veltype);
3352 0 : MDoppler mdop(velStart,DopType);
3353 0 : Quantity qrestfreq = restFreq.nelements() >0 ? restFreq[0]: Quantity(0.0, "Hz");
3354 0 : inStartFreq = MFrequency::fromDoppler(mdop, qrestfreq.getValue(Unit("Hz")), mfreqframe).getValue();
3355 : }
3356 : }
3357 : }
3358 :
3359 0 : return inStartFreq;
3360 :
3361 : }
3362 :
3363 0 : String SynthesisParamsImage::findSpecMode(const String& mode) const
3364 : {
3365 0 : String specmode;
3366 0 : specmode="channel";
3367 0 : 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 0 : if ( !(velStart.getValue()==0.0 && velStart.getUnit()=="" ) ||
3373 0 : !( velStep.getValue()==0.0)) {
3374 0 : specmode="velocity";
3375 : }
3376 0 : else if ( !(freqStart.getValue()==0.0 && freqStart.getUnit()=="") ||
3377 0 : !(freqStep.getValue()==0.0)) {
3378 0 : specmode="frequency";
3379 : }
3380 : }
3381 0 : return specmode;
3382 : }
3383 :
3384 :
3385 0 : Vector<Int> SynthesisParamsImage::decideNPolPlanes(const String& stokes) const
3386 : {
3387 0 : Vector<Int> whichStokes(0);
3388 0 : if(stokes=="I" || stokes=="Q" || stokes=="U" || stokes=="V" ||
3389 0 : stokes=="RR" ||stokes=="LL" ||
3390 0 : stokes=="XX" || stokes=="YY" ) {
3391 0 : whichStokes.resize(1);
3392 0 : whichStokes(0)=Stokes::type(stokes);
3393 : }
3394 0 : else if(stokes=="IV" || stokes=="IQ" ||
3395 0 : stokes=="RRLL" || stokes=="XXYY" ||
3396 0 : stokes=="QU" || stokes=="UV"){
3397 0 : whichStokes.resize(2);
3398 :
3399 0 : if(stokes=="IV"){ whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::V;}
3400 0 : else if(stokes=="IQ"){whichStokes[0]=Stokes::I; whichStokes[1]=Stokes::Q;}
3401 0 : else if(stokes=="RRLL"){whichStokes[0]=Stokes::RR; whichStokes[1]=Stokes::LL;}
3402 0 : else if(stokes=="XXYY"){whichStokes[0]=Stokes::XX; whichStokes[1]=Stokes::YY; }
3403 0 : 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 0 : 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 0 : else if(stokes=="IQUV"){
3416 0 : whichStokes.resize(4);
3417 0 : whichStokes(0)=Stokes::I; whichStokes(1)=Stokes::Q;
3418 0 : whichStokes(2)=Stokes::U; whichStokes(3)=Stokes::V;
3419 : }
3420 :
3421 0 : return whichStokes;
3422 : }// decidenpolplanes
3423 :
3424 0 : IPosition SynthesisParamsImage::shp() const
3425 : {
3426 0 : uInt nStokes = ( decideNPolPlanes(stokes) ).nelements();
3427 :
3428 0 : if( imsize[0]<=0 || imsize[1]<=0 || nStokes<=0 || nchan<=0 )
3429 : {
3430 0 : 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 0 : return IPosition( 4, imsize[0], imsize[1], nStokes, nchan );
3434 : }
3435 :
3436 0 : Record SynthesisParamsImage::getcsys() const
3437 : {
3438 0 : 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 0 : SynthesisParamsGrid::SynthesisParamsGrid():SynthesisParams()
3495 : {
3496 0 : setDefaults();
3497 0 : }
3498 :
3499 0 : SynthesisParamsGrid::~SynthesisParamsGrid()
3500 : {
3501 0 : }
3502 :
3503 :
3504 0 : void SynthesisParamsGrid::fromRecord(const Record &inrec)
3505 : {
3506 0 : setDefaults();
3507 :
3508 0 : String err("");
3509 :
3510 : try
3511 : {
3512 0 : err += readVal( inrec, String("imagename"), imageName);
3513 :
3514 : // FTMachine parameters
3515 0 : err += readVal( inrec, String("gridder"), gridder );
3516 0 : err += readVal( inrec, String("padding"), padding );
3517 0 : err += readVal( inrec, String("useautocorr"), useAutoCorr );
3518 0 : err += readVal( inrec, String("usedoubleprec"), useDoublePrec );
3519 0 : err += readVal( inrec, String("wprojplanes"), wprojplanes );
3520 0 : err += readVal( inrec, String("convfunc"), convFunc );
3521 :
3522 0 : err += readVal( inrec, String("vptable"), vpTable );
3523 :
3524 : //// convert 'gridder' to 'ftmachine' and 'mtype'
3525 0 : ftmachine="gridft";
3526 0 : mType="default";
3527 0 : if(gridder=="ft" || gridder=="gridft" || gridder=="standard" )
3528 0 : { ftmachine="gridft"; }
3529 0 : if( (gridder=="widefield" || gridder=="wproject" || gridder=="wprojectft" ) && (wprojplanes>1 || wprojplanes==-1))
3530 0 : { ftmachine="wprojectft";}
3531 :
3532 0 : if(gridder=="ftmosaic" || gridder=="mosaicft" || gridder=="mosaic" )
3533 0 : { ftmachine="mosaicft"; }
3534 0 : if(gridder=="imagemosaic") {
3535 0 : mType="imagemosaic";
3536 0 : if (wprojplanes>1 || wprojplanes==-1){ ftmachine="wprojectft"; }
3537 : }
3538 0 : if(gridder=="awproject" || gridder=="awprojectft" || gridder=="awp")
3539 0 : {ftmachine="awprojectft";}
3540 0 : if(gridder=="singledish") {
3541 0 : ftmachine="sd";
3542 : }
3543 :
3544 0 : String deconvolver;
3545 0 : err += readVal( inrec, String("deconvolver"), deconvolver );
3546 0 : if( deconvolver== "mtmfs" )
3547 0 : { mType="multiterm"; }// Takes precedence over imagemosaic
3548 :
3549 : // facets
3550 0 : err += readVal( inrec, String("facets"), facets);
3551 : // chanchunks
3552 0 : err += readVal( inrec, String("chanchunks"), chanchunks);
3553 :
3554 : // Spectral interpolation
3555 0 : err += readVal( inrec, String("interpolation"), interpolation );// not used in SI yet...
3556 : // Track moving source ?
3557 0 : err += readVal( inrec, String("distance"), distance );
3558 0 : err += readVal( inrec, String("tracksource"), trackSource );
3559 0 : err += readVal( inrec, String("trackdir"), trackDir );
3560 :
3561 : // The extra params for WB-AWP
3562 0 : err += readVal( inrec, String("aterm"), aTermOn );
3563 0 : err += readVal( inrec, String("psterm"), psTermOn );
3564 0 : err += readVal( inrec, String("mterm"), mTermOn );
3565 0 : err += readVal( inrec, String("wbawp"), wbAWP );
3566 0 : err += readVal( inrec, String("cfcache"), cfCache );
3567 0 : err += readVal( inrec, String("usepointing"), usePointing );
3568 0 : err += readVal( inrec, String("pointingoffsetsigdev"), pointingOffsetSigDev );
3569 0 : err += readVal( inrec, String("dopbcorr"), doPBCorr );
3570 0 : err += readVal( inrec, String("conjbeams"), conjBeams );
3571 0 : err += readVal( inrec, String("computepastep"), computePAStep );
3572 0 : err += readVal( inrec, String("rotatepastep"), rotatePAStep );
3573 :
3574 : // The extra params for single-dish
3575 0 : err += readVal( inrec, String("pointingcolumntouse"), pointingDirCol );
3576 0 : err += readVal( inrec, String("skypolthreshold"), skyPosThreshold );
3577 0 : err += readVal( inrec, String("convsupport"), convSupport );
3578 0 : err += readVal( inrec, String("truncate"), truncateSize );
3579 0 : err += readVal( inrec, String("gwidth"), gwidth );
3580 0 : err += readVal( inrec, String("jwidth"), jwidth );
3581 0 : err += readVal( inrec, String("minweight"), minWeight );
3582 0 : 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 0 : if( ftmachine=="awprojectft" && cfCache=="" )
3588 0 : {cfCache=imageName+".cf"; }
3589 :
3590 0 : if( ftmachine=="awprojectft" &&
3591 0 : usePointing==True &&
3592 0 : pointingOffsetSigDev.nelements() != 2 )
3593 : {
3594 : // Set the default to a large value so that it behaves like CASA 5.6's usepointing=True.
3595 0 : pointingOffsetSigDev.resize(2);
3596 0 : pointingOffsetSigDev[0]=600.0;
3597 0 : pointingOffsetSigDev[1]=600.0;
3598 : }
3599 :
3600 0 : err += verify();
3601 :
3602 : }
3603 0 : catch(AipsError &x)
3604 : {
3605 0 : err = err + x.getMesg() + "\n";
3606 : }
3607 :
3608 0 : if( err.length()>0 ) throw(AipsError("Invalid Gridding/FTM Parameter set : " + err));
3609 :
3610 0 : }
3611 :
3612 0 : String SynthesisParamsGrid::verify() const
3613 : {
3614 0 : String err;
3615 :
3616 : // Check for valid FTMachine type.
3617 : // Valid other params per FTM type, etc... ( check about nterms>1 )
3618 :
3619 0 : if( imageName=="" ) {err += "Please supply an image name\n";}
3620 :
3621 0 : if( (ftmachine != "gridft") && (ftmachine != "wprojectft") &&
3622 0 : (ftmachine != "mosaicft") && (ftmachine != "awprojectft") &&
3623 0 : (ftmachine != "mawprojectft") && (ftmachine != "protoft") &&
3624 0 : (ftmachine != "sd"))
3625 0 : { err += "Invalid ftmachine name. Must be one of 'gridft', 'wprojectft', 'mosaicft', 'awprojectft', 'mawpojectft'"; }
3626 :
3627 0 : if( ((ftmachine=="mosaicft") && (mType=="imagemosaic")) ||
3628 0 : ((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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : 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 0 : return err;
3667 : }
3668 :
3669 0 : void SynthesisParamsGrid::setDefaults()
3670 : {
3671 0 : imageName="";
3672 : // FTMachine parameters
3673 : //ftmachine="GridFT";
3674 0 : ftmachine="gridft";
3675 0 : gridder=ftmachine;
3676 0 : padding=1.2;
3677 0 : useAutoCorr=false;
3678 0 : useDoublePrec=true;
3679 0 : wprojplanes=1;
3680 0 : convFunc="SF";
3681 0 : vpTable="";
3682 :
3683 : // facets
3684 0 : facets=1;
3685 :
3686 : // chanchunks
3687 0 : chanchunks=1;
3688 :
3689 : // Spectral Axis interpolation
3690 0 : interpolation=String("nearest");
3691 :
3692 : //mosaic use pointing
3693 0 : usePointing=false;
3694 : // Moving phase center ?
3695 0 : distance=Quantity(0,"m");
3696 0 : trackSource=false;
3697 0 : trackDir=MDirection(Quantity(0.0, "deg"), Quantity(90.0, "deg"));
3698 :
3699 : // The extra params for WB-AWP
3700 0 : aTermOn = true;
3701 0 : psTermOn = true;
3702 0 : mTermOn = false;
3703 0 : wbAWP = true;
3704 0 : cfCache = "";
3705 0 : usePointing = false;
3706 0 : pointingOffsetSigDev.resize(0);
3707 : // pointingOffsetSigDev.set(30.0);
3708 0 : doPBCorr = true;
3709 0 : conjBeams = true;
3710 0 : computePAStep=360.0;
3711 0 : rotatePAStep=5.0;
3712 :
3713 : // extra params for single-dish
3714 0 : pointingDirCol = "";
3715 0 : skyPosThreshold = 0.0;
3716 0 : convSupport = -1;
3717 0 : truncateSize = Quantity(-1.0);
3718 0 : gwidth = Quantity(-1.0);
3719 0 : jwidth = Quantity(-1.0);
3720 0 : minWeight = 0.0;
3721 0 : clipMinMax = False;
3722 :
3723 : // Mapper type
3724 0 : mType = String("default");
3725 :
3726 0 : }
3727 :
3728 0 : Record SynthesisParamsGrid::toRecord() const
3729 : {
3730 0 : Record gridpar;
3731 :
3732 0 : gridpar.define("imagename", imageName);
3733 : // FTMachine params
3734 0 : gridpar.define("padding", padding);
3735 0 : gridpar.define("useautocorr",useAutoCorr );
3736 0 : gridpar.define("usedoubleprec", useDoublePrec);
3737 0 : gridpar.define("wprojplanes", wprojplanes);
3738 0 : gridpar.define("convfunc", convFunc);
3739 0 : gridpar.define("vptable", vpTable);
3740 :
3741 0 : gridpar.define("facets", facets);
3742 0 : gridpar.define("chanchunks", chanchunks);
3743 :
3744 0 : gridpar.define("interpolation",interpolation);
3745 :
3746 0 : gridpar.define("distance", QuantityToString(distance));
3747 0 : gridpar.define("tracksource", trackSource);
3748 0 : gridpar.define("trackdir", MDirectionToString( trackDir ));
3749 :
3750 0 : gridpar.define("aterm",aTermOn );
3751 0 : gridpar.define("psterm",psTermOn );
3752 0 : gridpar.define("mterm",mTermOn );
3753 0 : gridpar.define("wbawp", wbAWP);
3754 0 : gridpar.define("cfcache", cfCache);
3755 0 : gridpar.define("usepointing",usePointing );
3756 0 : gridpar.define("pointingoffsetsigdev", pointingOffsetSigDev);
3757 0 : gridpar.define("dopbcorr", doPBCorr);
3758 0 : gridpar.define("conjbeams",conjBeams );
3759 0 : gridpar.define("computepastep", computePAStep);
3760 0 : gridpar.define("rotatepastep", rotatePAStep);
3761 :
3762 0 : gridpar.define("pointingcolumntouse", pointingDirCol );
3763 0 : gridpar.define("skyposthreshold", skyPosThreshold );
3764 0 : gridpar.define("convsupport", convSupport );
3765 0 : gridpar.define("truncate", QuantityToString(truncateSize) );
3766 0 : gridpar.define("gwidth", QuantityToString(gwidth) );
3767 0 : gridpar.define("jwidth", QuantityToString(jwidth) );
3768 0 : gridpar.define("minweight", minWeight );
3769 0 : gridpar.define("clipminmax", clipMinMax );
3770 :
3771 0 : if( mType=="multiterm") gridpar.define("deconvolver","mtmfs");
3772 : /// else gridpar.define("deconvolver","singleterm");
3773 :
3774 0 : if( mType=="imagemosaic") gridpar.define("gridder","imagemosaic");
3775 0 : else gridpar.define("gridder", gridder);
3776 :
3777 : // gridpar.define("mtype", mType);
3778 :
3779 0 : return gridpar;
3780 : }
3781 :
3782 :
3783 :
3784 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////
3785 :
3786 : /////////////////////// Deconvolver Parameters
3787 :
3788 0 : SynthesisParamsDeconv::SynthesisParamsDeconv():SynthesisParams()
3789 : {
3790 0 : setDefaults();
3791 0 : }
3792 :
3793 0 : SynthesisParamsDeconv::~SynthesisParamsDeconv()
3794 : {
3795 0 : }
3796 :
3797 :
3798 0 : void SynthesisParamsDeconv::fromRecord(const Record &inrec)
3799 : {
3800 0 : setDefaults();
3801 :
3802 0 : String err("");
3803 :
3804 : try
3805 : {
3806 :
3807 0 : err += readVal( inrec, String("imagename"), imageName );
3808 0 : 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 0 : if( inrec.isDefined("startmodel") )
3814 : {
3815 0 : if( inrec.dataType("startmodel")==TpString )
3816 : {
3817 0 : String onemodel;
3818 0 : err += readVal( inrec, String("startmodel"), onemodel );
3819 0 : if( onemodel.length()>0 )
3820 : {
3821 0 : startModel.resize(1);
3822 0 : startModel[0] = onemodel;
3823 : }
3824 0 : else {startModel.resize();}
3825 : }
3826 0 : else if( inrec.dataType("startmodel")==TpArrayString ||
3827 0 : inrec.dataType("startmodel")==TpArrayBool)
3828 : {
3829 0 : 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 0 : err += readVal( inrec, String("id"), deconvolverId );
3838 0 : err += readVal( inrec, String("nterms"), nTaylorTerms );
3839 :
3840 0 : err += readVal( inrec, String("scales"), scales );
3841 0 : err += readVal( inrec, String("scalebias"), scalebias );
3842 :
3843 0 : err += readVal( inrec, String("usemask"), maskType );
3844 0 : if( maskType=="auto-thresh" )
3845 : {
3846 0 : autoMaskAlgorithm = "thresh";
3847 : }
3848 0 : else if( maskType=="auto-thesh2" )
3849 : {
3850 0 : autoMaskAlgorithm = "thresh2";
3851 : }
3852 0 : else if( maskType=="auto-onebox" )
3853 : {
3854 0 : autoMaskAlgorithm = "onebox";
3855 : }
3856 0 : else if( maskType=="user" || maskType=="pb" )
3857 : {
3858 0 : autoMaskAlgorithm = "";
3859 : }
3860 :
3861 :
3862 0 : if( inrec.isDefined("mask") )
3863 : {
3864 0 : if( inrec.dataType("mask")==TpString )
3865 : {
3866 0 : err+= readVal( inrec, String("mask"), maskString );
3867 : }
3868 0 : else if( inrec.dataType("mask")==TpArrayString )
3869 : {
3870 0 : err+= readVal( inrec, String("mask"), maskList );
3871 : }
3872 : }
3873 :
3874 0 : if( inrec.isDefined("pbmask") )
3875 : {
3876 0 : err += readVal( inrec, String("pbmask"), pbMask );
3877 : }
3878 0 : if( inrec.isDefined("maskthreshold") )
3879 : {
3880 0 : if( inrec.dataType("maskthreshold")==TpString )
3881 : {
3882 0 : err += readVal( inrec, String("maskthreshold"), maskThreshold );
3883 : //deal with the case a string is a float value without unit
3884 0 : Quantity testThresholdString;
3885 0 : Quantity::read(testThresholdString,maskThreshold);
3886 0 : if( testThresholdString.getUnit()=="" )
3887 : {
3888 0 : if(testThresholdString.getValue()<1.0)
3889 : {
3890 0 : fracOfPeak = testThresholdString.getValue();
3891 0 : 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 0 : if( inrec.isDefined("maskresolution") )
3913 : {
3914 0 : if( inrec.dataType("maskresolution")==TpString )
3915 : {
3916 0 : err += readVal(inrec, String("maskresolution"), maskResolution );
3917 : //deal with the case a string is a float value without unit
3918 0 : Quantity testResolutionString;
3919 0 : Quantity::read(testResolutionString,maskResolution);
3920 0 : if( testResolutionString.getUnit()=="" )
3921 : {
3922 0 : maskResByBeam = testResolutionString.getValue();
3923 0 : 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 0 : if( inrec.isDefined("nmask") )
3939 : {
3940 0 : if( inrec.dataType("nmask")==TpInt )
3941 : {
3942 0 : err+= readVal(inrec, String("nmask"), nMask );
3943 : }
3944 : else
3945 : {
3946 0 : err+= "nmask must be an integer\n";
3947 : }
3948 : }
3949 0 : if( inrec.isDefined("autoadjust") )
3950 : {
3951 0 : if( inrec.dataType("autoadjust")==TpBool )
3952 : {
3953 0 : 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 0 : if (inrec.isDefined("fusedthreshold"))
3962 : {
3963 0 : if (inrec.dataType("fusedthreshold") == TpFloat || inrec.dataType("fusedthreshold") == TpDouble)
3964 0 : err += readVal(inrec, String("fusedthreshold"), fusedThreshold);
3965 : else
3966 0 : err += "fusedthreshold must be a float or double";
3967 : }
3968 0 : if (inrec.isDefined("specmode"))
3969 : {
3970 0 : if(inrec.dataType("specmode") == TpString)
3971 0 : 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 0 : if (inrec.isDefined("largestscale"))
3977 : {
3978 0 : if (inrec.dataType("largestscale") == TpInt)
3979 0 : 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 0 : if( inrec.isDefined("sidelobethreshold"))
3985 : {
3986 0 : if(inrec.dataType("sidelobethreshold")==TpFloat || inrec.dataType("sidelobethreshold")==TpDouble )
3987 : {
3988 0 : err+= readVal(inrec, String("sidelobethreshold"), sidelobeThreshold );
3989 : }
3990 : else
3991 : {
3992 0 : err+= "sidelobethreshold must be a float or double";
3993 : }
3994 : }
3995 :
3996 0 : if( inrec.isDefined("noisethreshold"))
3997 : {
3998 0 : if(inrec.dataType("noisethreshold")==TpFloat || inrec.dataType("noisethreshold")==TpDouble )
3999 : {
4000 0 : err+= readVal(inrec, String("noisethreshold"), noiseThreshold );
4001 : }
4002 : else
4003 : {
4004 0 : err+= "noisethreshold must be a float or double";
4005 : }
4006 : }
4007 0 : if( inrec.isDefined("lownoisethreshold"))
4008 : {
4009 0 : if(inrec.dataType("lownoisethreshold")==TpFloat || inrec.dataType("lownoisethreshold")==TpDouble )
4010 : {
4011 0 : err+= readVal(inrec, String("lownoisethreshold"), lowNoiseThreshold );
4012 : }
4013 : else
4014 : {
4015 0 : err+= "lownoisethreshold must be a float or double";
4016 : }
4017 : }
4018 0 : if( inrec.isDefined("negativethreshold"))
4019 : {
4020 0 : if(inrec.dataType("negativethreshold")==TpFloat || inrec.dataType("negativethreshold")==TpDouble )
4021 : {
4022 0 : err+= readVal(inrec, String("negativethreshold"), negativeThreshold );
4023 : }
4024 : else
4025 : {
4026 0 : err+= "negativethreshold must be a float or double";
4027 : }
4028 : }
4029 0 : if( inrec.isDefined("smoothfactor"))
4030 : {
4031 0 : if( inrec.dataType("smoothfactor")==TpFloat || inrec.dataType("smoothfactor")==TpDouble )
4032 : {
4033 0 : err+= readVal(inrec, String("smoothfactor"), smoothFactor );
4034 : }
4035 : else
4036 : {
4037 0 : err+= "smoothfactor must be a float or double";
4038 : }
4039 : }
4040 0 : if( inrec.isDefined("minbeamfrac"))
4041 : {
4042 0 : if( inrec.dataType("minbeamfrac")==TpFloat || inrec.dataType("minbeamfrac")==TpDouble )
4043 : {
4044 0 : 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 0 : if( inrec.isDefined("cutthreshold"))
4058 : {
4059 0 : if( inrec.dataType("cutthreshold")==TpFloat || inrec.dataType("cutthreshold")==TpDouble )
4060 : {
4061 0 : err+= readVal(inrec, String("cutthreshold"), cutThreshold );
4062 : }
4063 : else {
4064 0 : err+= "cutthreshold must be a float or double";
4065 : }
4066 : }
4067 0 : if( inrec.isDefined("growiterations"))
4068 : {
4069 0 : if (inrec.dataType("growiterations")==TpInt) {
4070 0 : err+= readVal(inrec, String("growiterations"), growIterations );
4071 : }
4072 : else {
4073 0 : err+= "growiterations must be an integer\n";
4074 : }
4075 : }
4076 0 : if( inrec.isDefined("dogrowprune"))
4077 : {
4078 0 : if (inrec.dataType("dogrowprune")==TpBool) {
4079 0 : err+= readVal(inrec, String("dogrowprune"), doGrowPrune );
4080 : }
4081 : else {
4082 0 : err+= "dogrowprune must be a bool\n";
4083 : }
4084 : }
4085 0 : if( inrec.isDefined("minpercentchange"))
4086 : {
4087 0 : if (inrec.dataType("minpercentchange")==TpFloat || inrec.dataType("minpercentchange")==TpDouble ) {
4088 0 : err+= readVal(inrec, String("minpercentchange"), minPercentChange );
4089 : }
4090 : else {
4091 0 : err+= "minpercentchange must be a float or double";
4092 : }
4093 : }
4094 0 : if( inrec.isDefined("verbose"))
4095 : {
4096 0 : if (inrec.dataType("verbose")==TpBool ) {
4097 0 : err+= readVal(inrec, String("verbose"), verbose);
4098 : }
4099 : else {
4100 0 : err+= "verbose must be a bool";
4101 : }
4102 : }
4103 0 : if( inrec.isDefined("fastnoise"))
4104 : {
4105 0 : if (inrec.dataType("fastnoise")==TpBool ) {
4106 0 : err+= readVal(inrec, String("fastnoise"), fastnoise);
4107 : }
4108 : else {
4109 0 : err+= "fastnoise must be a bool";
4110 : }
4111 : }
4112 0 : if( inrec.isDefined("nsigma") )
4113 : {
4114 0 : if(inrec.dataType("nsigma")==TpFloat || inrec.dataType("nsigma")==TpDouble ) {
4115 0 : 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 0 : if( inrec.isDefined("noRequireSumwt") )
4128 : {
4129 0 : if (inrec.dataType("noRequireSumwt")==TpBool) {
4130 0 : err+= readVal(inrec, String("noRequireSumwt"), noRequireSumwt);
4131 : }
4132 : else {
4133 0 : err+= "noRequireSumwt must be a bool";
4134 : }
4135 : }
4136 0 : if( inrec.isDefined("fullsummary") )
4137 : {
4138 0 : if (inrec.dataType("fullsummary")==TpBool) {
4139 0 : err+= readVal(inrec, String("fullsummary"), fullsummary);
4140 : }
4141 : else {
4142 0 : err+= "fullsummary must be a bool";
4143 : }
4144 : }
4145 0 : if( inrec.isDefined("restoringbeam") )
4146 : {
4147 0 : String errinfo("");
4148 : try {
4149 :
4150 0 : if( inrec.dataType("restoringbeam")==TpString )
4151 : {
4152 0 : 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 0 : if( (! usebeam.matches("common")) && usebeam.length()!=0 )
4157 : {
4158 0 : Quantity bsize;
4159 0 : err += readVal( inrec, String("restoringbeam"), bsize );
4160 0 : restoringbeam.setMajorMinor( bsize, bsize );
4161 0 : usebeam = String("");
4162 : }
4163 0 : errinfo = usebeam;
4164 : }
4165 0 : 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 0 : if( inrec.isDefined("interactive") )
4201 : {
4202 0 : if( inrec.dataType("interactive")==TpBool )
4203 0 : {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 0 : err += verify();
4211 :
4212 : }
4213 0 : catch(AipsError &x)
4214 : {
4215 0 : err = err + x.getMesg() + "\n";
4216 : }
4217 :
4218 0 : if( err.length()>0 ) throw(AipsError("Invalid Deconvolver Parameter set : " + err));
4219 :
4220 0 : }
4221 :
4222 0 : String SynthesisParamsDeconv::verify() const
4223 : {
4224 0 : String err;
4225 :
4226 0 : 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 0 : if( maskString.length()>0 )
4230 : {
4231 0 : File fp( imageName+".mask" );
4232 0 : if( fp.exists() ) err += "Mask image " + imageName+".mask exists, but a specific input mask of " + maskString + " has also been supplied. Please either reset mask='' to reuse the existing mask, or delete " + imageName + ".mask before restarting";
4233 : }
4234 :
4235 0 : if( pbMask >= 1.0)
4236 0 : {err += "pbmask must be < 1.0 \n"; }
4237 0 : else if( pbMask < 0.0)
4238 0 : {err += "pbmask must be a positive value \n"; }
4239 :
4240 0 : 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 0 : if ( fracOfPeak >= 1.0)
4249 0 : {err += "fracofpeak must be < 1.0 \n"; }
4250 0 : else if ( fracOfPeak < 0.0)
4251 0 : {err += "fracofpeak must be a positive value \n"; }
4252 :
4253 0 : return err;
4254 : }
4255 :
4256 0 : void SynthesisParamsDeconv::setDefaults()
4257 : {
4258 0 : imageName="";
4259 0 : algorithm="hogbom";
4260 0 : startModel=Vector<String>(0);
4261 0 : deconvolverId=0;
4262 0 : nTaylorTerms=1;
4263 0 : scales.resize(1); scales[0]=0.0;
4264 0 : scalebias=0.6;
4265 0 : maskType="none";
4266 0 : maskString="";
4267 0 : maskList.resize(1); maskList[0]="";
4268 0 : pbMask=0.0;
4269 0 : autoMaskAlgorithm="thresh";
4270 0 : maskThreshold="";
4271 0 : maskResolution="";
4272 0 : fracOfPeak=0.0;
4273 0 : nMask=0;
4274 0 : interactive=false;
4275 0 : autoAdjust=False;
4276 0 : fusedThreshold = 0.0;
4277 0 : specmode="mfs";
4278 0 : largestscale = -1;
4279 0 : }
4280 :
4281 0 : Record SynthesisParamsDeconv::toRecord() const
4282 : {
4283 0 : Record decpar;
4284 :
4285 0 : decpar.define("imagename", imageName);
4286 0 : decpar.define("deconvolver", algorithm);
4287 0 : decpar.define("startmodel",startModel);
4288 0 : decpar.define("id",deconvolverId);
4289 0 : decpar.define("nterms",nTaylorTerms);
4290 0 : decpar.define("scales",scales);
4291 0 : decpar.define("scalebias",scalebias);
4292 0 : decpar.define("usemask",maskType);
4293 0 : decpar.define("fusedthreshold", fusedThreshold);
4294 0 : decpar.define("specmode", specmode);
4295 0 : decpar.define("largestscale", largestscale);
4296 0 : if( maskList.nelements()==1 && maskList[0]=="")
4297 : {
4298 0 : decpar.define("mask",maskString);
4299 : }
4300 : else {
4301 0 : decpar.define("mask",maskList);
4302 : }
4303 0 : decpar.define("pbmask",pbMask);
4304 0 : if (fracOfPeak > 0.0)
4305 : {
4306 0 : decpar.define("maskthreshold",fracOfPeak);
4307 : }
4308 : else
4309 : {
4310 0 : decpar.define("maskthreshold",maskThreshold);
4311 : }
4312 0 : decpar.define("maskresolution",maskResolution);
4313 0 : decpar.define("nmask",nMask);
4314 0 : decpar.define("autoadjust",autoAdjust);
4315 0 : decpar.define("sidelobethreshold",sidelobeThreshold);
4316 0 : decpar.define("noisethreshold",noiseThreshold);
4317 0 : decpar.define("lownoisethreshold",lowNoiseThreshold);
4318 0 : decpar.define("negativethreshold",negativeThreshold);
4319 0 : decpar.define("smoothfactor",smoothFactor);
4320 0 : decpar.define("minbeamfrac",minBeamFrac);
4321 0 : decpar.define("cutthreshold",cutThreshold);
4322 0 : decpar.define("growiterations",growIterations);
4323 0 : decpar.define("dogrowprune",doGrowPrune);
4324 0 : decpar.define("minpercentchange",minPercentChange);
4325 0 : decpar.define("verbose", verbose);
4326 0 : decpar.define("fastnoise", fastnoise);
4327 0 : decpar.define("interactive",interactive);
4328 0 : decpar.define("nsigma",nsigma);
4329 0 : decpar.define("noRequireSumwt",noRequireSumwt);
4330 0 : decpar.define("fullsummary",fullsummary);
4331 :
4332 0 : return decpar;
4333 : }
4334 :
4335 : /////////////////////////////////////////////////////////////////////////////////////////////////////
4336 :
4337 :
4338 : } //# NAMESPACE CASA - END
4339 :
|
