Line data Source code
1 : //# newsimulator.cc: Simulation program
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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: Simulator.cc,v 1.1.2.4 2006/10/06 21:03:19 kgolap Exp $
27 :
28 : #include <stdexcept>
29 : #include <casacore/casa/Arrays/Matrix.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <iostream>
33 :
34 : #include <casacore/casa/Logging.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/casa/Containers/Record.h>
38 : #include <casacore/casa/Containers/RecordInterface.h>
39 :
40 : #include <casacore/tables/TaQL/TableParse.h>
41 : #include <casacore/tables/Tables/TableRecord.h>
42 : #include <casacore/tables/Tables/TableDesc.h>
43 : #include <casacore/tables/Tables/TableLock.h>
44 : #include <casacore/tables/TaQL/ExprNode.h>
45 :
46 : #include <casacore/casa/BasicSL/String.h>
47 : #include <casacore/casa/Utilities/Assert.h>
48 : #include <casacore/casa/Utilities/Fallible.h>
49 :
50 : #include <casacore/casa/BasicSL/Constants.h>
51 :
52 : #include <casacore/casa/Logging/LogSink.h>
53 : #include <casacore/casa/Logging/LogMessage.h>
54 :
55 : #include <casacore/casa/Arrays/ArrayMath.h>
56 :
57 : #include <msvis/MSVis/VisSet.h>
58 : #include <msvis/MSVis/VisSetUtil.h>
59 : #include <synthesis/MeasurementComponents/VisCal.h>
60 : #include <synthesis/MeasurementComponents/VisCalGlobals.h>
61 : #include <casacore/ms/MSOper/NewMSSimulator.h>
62 :
63 : #include <casacore/measures/Measures/Stokes.h>
64 : #include <casacore/casa/Quanta/UnitMap.h>
65 : #include <casacore/casa/Quanta/UnitVal.h>
66 : #include <casacore/casa/Quanta/MVAngle.h>
67 : #include <casacore/measures/Measures/MDirection.h>
68 : #include <casacore/measures/Measures/MPosition.h>
69 : #include <casacore/casa/Quanta/MVEpoch.h>
70 : #include <casacore/measures/Measures/MEpoch.h>
71 :
72 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
73 : #include <casacore/ms/MeasurementSets/MSColumns.h>
74 :
75 : #include <casacore/ms/MSOper/MSSummary.h>
76 : #include <synthesis/MeasurementEquations/SkyEquation.h>
77 : #include <synthesis/MeasurementComponents/ImageSkyModel.h>
78 : #include <synthesis/MeasurementComponents/SimACohCalc.h>
79 : #include <synthesis/MeasurementComponents/SimACoh.h>
80 : //#include <synthesis/MeasurementComponents/SimVisJones.h>
81 : #include <synthesis/TransformMachines/VPSkyJones.h>
82 : #include <synthesis/TransformMachines/StokesImageUtil.h>
83 : #include <casacore/lattices/LEL/LatticeExpr.h>
84 :
85 : #include <synthesis/MeasurementEquations/Simulator.h>
86 : #include <synthesis/MeasurementComponents/CleanImageSkyModel.h>
87 : #include <synthesis/MeasurementComponents/GridBoth.h>
88 : #include <synthesis/TransformMachines/WProjectFT.h>
89 : #include <synthesis/MeasurementComponents/GridBoth.h>
90 : #include <synthesis/TransformMachines/MosaicFTNew.h>
91 : #include <synthesis/MeasurementEquations/VPManager.h>
92 : #include <synthesis/TransformMachines/HetArrayConvFunc.h> //2016
93 : #include <synthesis/TransformMachines/SimpleComponentFTMachine.h>
94 : #include <casacore/casa/OS/HostInfo.h>
95 : #include <casacore/images/Images/PagedImage.h>
96 : #include <casacore/casa/Arrays/Cube.h>
97 : #include <casacore/casa/Arrays/Vector.h>
98 : #include <sstream>
99 :
100 : #include <casacore/casa/namespace.h>
101 :
102 : namespace casa {
103 :
104 0 : Simulator::Simulator():
105 : msname_p(String("")), ms_p(0), mssel_p(0), vs_p(0),
106 : seed_p(11111),
107 : ac_p(0), distance_p(0), ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
108 : sim_p(0),
109 : // epJ_p(0),
110 : epJTableName_p(),
111 0 : prtlev_(0)
112 : {
113 0 : }
114 :
115 :
116 29 : Simulator::Simulator(String& msname)
117 : : msname_p(msname), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
118 : ac_p(0), distance_p(0),ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
119 : sim_p(0),
120 : // epJ_p(0),
121 : epJTableName_p(),
122 29 : prtlev_(0)
123 : {
124 29 : defaults();
125 :
126 29 : if(!sim_p) {
127 29 : sim_p= new NewMSSimulator(msname);
128 : //sim_p->setPrtlev(prtlev());
129 : }
130 :
131 29 : ve_p.setPrtlev(prtlev());
132 :
133 : // Make a MeasurementSet object for the disk-base MeasurementSet that we just
134 : // created
135 29 : ms_p = new MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
136 29 : Table::Update);
137 29 : AlwaysAssert(ms_p, AipsError);
138 :
139 29 : }
140 :
141 :
142 16 : Simulator::Simulator(MeasurementSet &theMs)
143 : : msname_p(""), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
144 : ac_p(0), distance_p(0), ve_p(), vc_p(), vp_p(0), gvp_p(0),
145 : sim_p(0),
146 : // epJ_p(0),
147 : epJTableName_p(),
148 16 : prtlev_(0)
149 : {
150 48 : LogIO os(LogOrigin("simulator", "simulator(MeasurementSet& theMs)"));
151 :
152 16 : defaults();
153 :
154 16 : msname_p=theMs.tableName();
155 :
156 16 : if(!sim_p) {
157 16 : sim_p= new NewMSSimulator(theMs);
158 : //sim_p->setPrtlev(prtlev());
159 : }
160 :
161 16 : ve_p.setPrtlev(prtlev());
162 :
163 16 : ms_p = new MeasurementSet(theMs);
164 16 : AlwaysAssert(ms_p, AipsError);
165 :
166 : // get info from the MS into Simulator:
167 16 : if (!getconfig())
168 0 : os << "Can't find antenna information for loaded MS" << LogIO::WARN;
169 16 : if (!sim_p->getSpWindows(nSpw,spWindowName_p,nChan_p,startFreq_p,freqInc_p,stokesString_p))
170 0 : os << "Can't find spectral window information for loaded MS" << LogIO::WARN;
171 16 : freqRes_p.resize(nSpw);
172 32 : for (Int i=0;i<nSpw;++i)
173 16 : freqRes_p[i]=freqInc_p[i];
174 16 : if (!sim_p->getFields(nField,sourceName_p,sourceDirection_p,calCode_p))
175 0 : os << "Can't find Field/Source information for loaded MS" << LogIO::WARN;
176 :
177 16 : if (!sim_p->getFeedMode(feedMode_p))
178 0 : os << "Can't find Feed information for loaded MS" << LogIO::WARN;
179 : else
180 16 : feedsHaveBeenSet=true;
181 :
182 16 : }
183 :
184 :
185 :
186 0 : Simulator::Simulator(const Simulator &other)
187 : : msname_p(""), ms_p(0), vs_p(0), seed_p(11111),
188 : ac_p(0), distance_p(0),ve_p(), vc_p(), vp_p(0), gvp_p(0),
189 : sim_p(0),
190 : // epJ_p(0),
191 : epJTableName_p(),
192 0 : prtlev_(0)
193 : {
194 0 : defaults();
195 0 : ms_p = new MeasurementSet(*other.ms_p);
196 0 : if(other.mssel_p) {
197 0 : mssel_p = new MeasurementSet(*other.mssel_p);
198 : }
199 0 : }
200 :
201 0 : Simulator &Simulator::operator=(const Simulator &other)
202 : {
203 0 : if (ms_p && this != &other) {
204 0 : *ms_p = *(other.ms_p);
205 : }
206 0 : if (mssel_p && this != &other && other.mssel_p) {
207 0 : *mssel_p = *(other.mssel_p);
208 : }
209 0 : if (vs_p && this != &other) {
210 0 : *vs_p = *(other.vs_p);
211 : }
212 0 : if (ac_p && this != &other) {
213 0 : *ac_p = *(other.ac_p);
214 : }
215 :
216 : // TBD VisEquation/VisCal stuff
217 :
218 0 : if (vp_p && this != &other) {
219 0 : *vp_p = *(other.vp_p);
220 : }
221 0 : if (gvp_p && this != &other) {
222 0 : *gvp_p = *(other.gvp_p);
223 : }
224 0 : if (sim_p && this != &other) {
225 0 : *sim_p = *(other.sim_p);
226 : }
227 : // if (epJ_p && this != &other) *epJ_p = *(other.epJ_p);
228 0 : return *this;
229 : }
230 :
231 45 : Simulator::~Simulator()
232 : {
233 45 : if (ms_p) {
234 45 : ms_p->relinquishAutoLocks();
235 45 : ms_p->unlock();
236 45 : delete ms_p;
237 : }
238 45 : ms_p = 0;
239 45 : if (mssel_p) {
240 43 : mssel_p->relinquishAutoLocks();
241 43 : mssel_p->unlock();
242 43 : delete mssel_p;
243 : }
244 45 : mssel_p = 0;
245 45 : if (vs_p) {
246 43 : delete vs_p;
247 : }
248 45 : vs_p = 0;
249 :
250 : // Delete all vis-plane calibration corruption terms
251 45 : resetviscal();
252 :
253 : // Delete all im-plane calibration corruption terms
254 45 : resetimcal();
255 :
256 45 : if(sim_p) delete sim_p; sim_p = 0;
257 :
258 45 : if(sm_p) delete sm_p; sm_p = 0;
259 45 : if(ft_p) delete ft_p; ft_p = 0;
260 45 : if(cft_p) delete cft_p; cft_p = 0;
261 :
262 45 : }
263 :
264 :
265 45 : void Simulator::defaults()
266 : {
267 45 : UnitMap::putUser("Pixel", UnitVal(1.0), "Pixel solid angle");
268 45 : UnitMap::putUser("Beam", UnitVal(1.0), "Beam solid angle");
269 45 : gridfunction_p="SF";
270 : // Use half the machine memory as cache. The user can override
271 : // this via the setoptions function().
272 45 : cache_p=(HostInfo::memoryTotal()/8)*1024;
273 :
274 45 : tile_p=16;
275 45 : ftmachine_p="gridft";
276 45 : padding_p=1.3;
277 45 : facets_p=1;
278 45 : sm_p = 0;
279 45 : ft_p = 0;
280 45 : cft_p = 0;
281 45 : vp_p = 0;
282 45 : gvp_p = 0;
283 45 : sim_p = 0;
284 45 : images_p = 0;
285 45 : nmodels_p = 1;
286 : // info for configurations
287 45 : areStationCoordsSet_p = false;
288 45 : telescope_p = "UNSET";
289 45 : nmodels_p = 0;
290 :
291 : // info for fields and schedule:
292 45 : nField=0;
293 45 : sourceName_p.resize(1);
294 45 : sourceName_p[0]="UNSET";
295 45 : calCode_p.resize(1);
296 45 : calCode_p[0]="";
297 45 : sourceDirection_p.resize(1);
298 :
299 : // info for spectral windows
300 45 : nSpw=0;
301 45 : spWindowName_p.resize(1);
302 45 : nChan_p.resize(1);
303 45 : startFreq_p.resize(1);
304 45 : freqInc_p.resize(1);
305 45 : freqRes_p.resize(1);
306 45 : stokesString_p.resize(1);
307 45 : spWindowName_p[0]="UNSET";
308 45 : nChan_p[0]=1;
309 45 : startFreq_p[0]=Quantity(50., "GHz");
310 45 : freqInc_p[0]=Quantity(0.1, "MHz");
311 45 : freqRes_p[0]=Quantity(0.1, "MHz");
312 45 : stokesString_p[0]="RR RL LR LL";
313 :
314 : // feeds
315 45 : feedMode_p = "perfect R L";
316 45 : nFeeds_p = 1;
317 45 : feedsHaveBeenSet = false;
318 45 : feedsInitialized = false;
319 :
320 : // times
321 45 : integrationTime_p = Quantity(10.0, "s");
322 45 : useHourAngle_p=true;
323 45 : refTime_p = MEpoch(Quantity(0.0, "s"), MEpoch::UTC);
324 45 : timesHaveBeenSet_p=false;
325 :
326 : // VP stuff
327 45 : doVP_p=false;
328 45 : doDefaultVP_p = true;
329 :
330 45 : };
331 :
332 :
333 0 : Bool Simulator::close()
334 : {
335 0 : LogIO os(LogOrigin("Simulator", "close()", WHERE));
336 : os << "Closing MeasurementSet and detaching from Simulator"
337 0 : << LogIO::POST;
338 :
339 : // Delete all im-plane calibration corruption terms
340 0 : resetimcal();
341 : // Delete all vis-plane calibration corruption terms
342 0 : resetviscal();
343 :
344 0 : ms_p->unlock();
345 0 : if(mssel_p) mssel_p->unlock();
346 0 : if(vs_p) delete vs_p; vs_p = 0;
347 0 : if(mssel_p) delete mssel_p; mssel_p = 0;
348 0 : if(ms_p) delete ms_p; ms_p = 0;
349 0 : if(sm_p) delete sm_p; sm_p = 0;
350 0 : if(ft_p) delete ft_p; ft_p = 0;
351 0 : if(cft_p) delete cft_p; cft_p = 0;
352 :
353 0 : return true;
354 : }
355 :
356 45 : Bool Simulator::resetviscal() {
357 135 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
358 : try {
359 :
360 45 : os << "Resetting all visibility corruption components" << LogIO::POST;
361 :
362 : // The noise term (for now)
363 45 : if(ac_p) delete ac_p; ac_p=0;
364 :
365 : // Delete all VisCals
366 77 : for (uInt i=0;i<vc_p.nelements();++i)
367 32 : if (vc_p[i]) delete vc_p[i];
368 45 : vc_p.resize(0,true);
369 :
370 : // reset the VisEquation (by sending an empty vc_p)
371 45 : ve_p.setapply(vc_p);
372 :
373 0 : } catch (AipsError x) {
374 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
375 0 : return false;
376 : }
377 45 : return true;
378 : }
379 :
380 :
381 45 : Bool Simulator::resetimcal() {
382 135 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
383 : try {
384 :
385 45 : os << "Reset all image-plane corruption components" << LogIO::POST;
386 :
387 45 : if(vp_p) delete vp_p; vp_p=0;
388 45 : if(gvp_p) delete gvp_p; gvp_p=0;
389 : /*
390 : // if(epJ_p) delete epJ_p; epJ_p=0;
391 : */
392 :
393 0 : } catch (AipsError x) {
394 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
395 0 : return false;
396 : }
397 45 : return true;
398 : }
399 :
400 :
401 0 : Bool Simulator::reset() {
402 0 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
403 : try {
404 :
405 : // reset vis-plane cal terms
406 0 : resetviscal();
407 :
408 : // reset im-plane cal terms
409 0 : resetimcal();
410 :
411 0 : } catch (AipsError x) {
412 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
413 0 : return false;
414 : }
415 0 : return true;
416 : }
417 :
418 :
419 0 : String Simulator::name() const
420 : {
421 0 : if (detached()) {
422 0 : return "none";
423 : }
424 0 : return msname_p;
425 : }
426 :
427 0 : String Simulator::state()
428 : {
429 0 : ostringstream os;
430 0 : os << "Need to write the state() method!" << LogIO::POST;
431 0 : if(doVP_p) {
432 0 : os << " Primary beam correction is enabled" << endl;
433 : }
434 0 : return String(os);
435 : }
436 :
437 0 : Bool Simulator::summary()
438 : {
439 0 : LogIO os(LogOrigin("Simulator", "summary()", WHERE));
440 0 : createSummary(os);
441 0 : predictSummary(os);
442 0 : corruptSummary(os);
443 :
444 0 : return true;
445 : }
446 :
447 :
448 0 : Bool Simulator::createSummary(LogIO& os)
449 : {
450 0 : Bool configResult = configSummary(os);
451 0 : Bool fieldResult = fieldSummary(os);
452 0 : Bool windowResult = spWindowSummary(os);
453 0 : Bool feedResult = feedSummary(os);
454 :
455 0 : if (!configResult && !fieldResult && !windowResult && !feedResult) {
456 0 : os << "=======================================" << LogIO::POST;
457 0 : os << "No create-type information has been set" << LogIO::POST;
458 0 : os << "=======================================" << LogIO::POST;
459 0 : return false;
460 : } else {
461 : // user has set at least ONE, so we report on each
462 0 : if (!configResult) {
463 0 : os << "No configuration information set yet, but other create-type info HAS been set" << LogIO::POST;
464 : }
465 0 : if (!fieldResult) {
466 0 : os << "No field information set yet, but other create-type info HAS been set" << LogIO::POST;
467 : }
468 0 : if (!windowResult) {
469 0 : os << "No window information set yet, but other create-type info HAS been set" << LogIO::POST;
470 : }
471 0 : if (!feedResult) {
472 0 : os << "No feed information set yet, but other create-type info HAS been set" << LogIO::POST;
473 0 : os << "(feeds will default to perfect R-L feeds if not set)" << LogIO::POST;
474 : }
475 0 : os << "======================================================================" << LogIO::POST;
476 : }
477 0 : return true;
478 : }
479 :
480 :
481 :
482 0 : Bool Simulator::configSummary(LogIO& os)
483 : {
484 0 : if ( ! areStationCoordsSet_p ) {
485 0 : return false;
486 : } else {
487 0 : os << "----------------------------------------------------------------------" << LogIO::POST;
488 0 : os << "Generating (u,v,w) using this configuration: " << LogIO::POST;
489 0 : os << " x y z diam mount station " << LogIO::POST;
490 0 : for (uInt i=0; i< x_p.nelements(); i++) {
491 0 : os << x_p(i)
492 0 : << " " << y_p(i)
493 0 : << " " << z_p(i)
494 0 : << " " << diam_p(i)
495 0 : << " " << mount_p(i)
496 0 : << " " << padName_p(i)
497 0 : << LogIO::POST;
498 : }
499 0 : os << " Coordsystem = " << coordsystem_p << LogIO::POST;
500 : os << " RefLocation = " <<
501 0 : mRefLocation_p.getAngle("deg").getValue("deg") << LogIO::POST;
502 : }
503 0 : return true;
504 :
505 : }
506 :
507 :
508 :
509 0 : Bool Simulator::fieldSummary(LogIO& os)
510 : {
511 0 : os << "----------------------------------------------------------------------" << LogIO::POST;
512 0 : os << " Field information: " << LogIO::POST;
513 0 : if (nField==0)
514 0 : os << "NO Field window information set" << LogIO::POST;
515 : else
516 0 : os << " Name direction calcode" << LogIO::POST;
517 0 : for (Int i=0;i<nField;i++)
518 0 : os << sourceName_p[i]
519 0 : << " " << formatDirection(sourceDirection_p[i])
520 0 : << " " << calCode_p[i]
521 0 : << LogIO::POST;
522 0 : return true;
523 : }
524 :
525 :
526 :
527 0 : Bool Simulator::timeSummary(LogIO& os)
528 : {
529 0 : if(integrationTime_p.getValue("s") <= 0.0) {
530 0 : return false;
531 : } else {
532 0 : os << "----------------------------------------------------------------------" << LogIO::POST;
533 0 : os << " Time information: " << LogIO::POST;
534 : os << " integration time = " << integrationTime_p.getValue("s")
535 0 : << " s" << LogIO::POST;
536 0 : os << " reference time = " << MVTime(refTime_p.get("s").getValue("d")).string()
537 0 : << LogIO::POST;
538 : }
539 0 : return true;
540 : }
541 :
542 :
543 :
544 0 : Bool Simulator::spWindowSummary(LogIO& os)
545 : {
546 0 : os << "----------------------------------------------------------------------" << LogIO::POST;
547 0 : os << " Spectral Windows information: " << LogIO::POST;
548 0 : if (nSpw==0)
549 0 : os << "NO Spectral window information set" << LogIO::POST;
550 : else
551 0 : os << " Name nchan freq[GHz] freqInc[MHz] freqRes[MHz] stokes" << LogIO::POST;
552 0 : for (Int i=0;i<nSpw;i++)
553 0 : os << spWindowName_p[i]
554 0 : << " " << nChan_p[i]
555 0 : << " " << startFreq_p[i].getValue("GHz")
556 0 : << " " << freqInc_p[i].getValue("MHz")
557 0 : << " " << freqRes_p[i].getValue("MHz")
558 0 : << " " << stokesString_p[i]
559 0 : << LogIO::POST;
560 0 : return true;
561 : }
562 :
563 :
564 0 : Bool Simulator::feedSummary(LogIO& os)
565 : {
566 0 : if (!feedsHaveBeenSet) {
567 0 : return false;
568 : } else {
569 0 : os << "----------------------------------------------------------------------" << LogIO::POST;
570 0 : os << " Feed information: " << LogIO::POST;
571 0 : os << feedMode_p << LogIO::POST;
572 : }
573 0 : return true;
574 : }
575 :
576 :
577 0 : Bool Simulator::predictSummary(LogIO& os)
578 : {
579 0 : Bool vpResult = vpSummary(os);
580 0 : Bool optionsResult = optionsSummary(os);
581 :
582 : // keep compiler happy
583 0 : if (!vpResult && !optionsResult) {}
584 0 : return true;
585 : }
586 :
587 :
588 0 : Bool Simulator::vpSummary(LogIO& /*os*/)
589 : {
590 0 : if (vp_p) {
591 0 : vp_p->summary();
592 0 : return true;
593 : } else {
594 0 : return false;
595 : }
596 : }
597 :
598 :
599 0 : Bool Simulator::optionsSummary(LogIO& /*os*/)
600 : {
601 0 : return true;
602 : }
603 :
604 :
605 0 : Bool Simulator::corruptSummary(LogIO& os)
606 : {
607 0 : if (vc_p.nelements()<1 && !ac_p) {
608 0 : os << "===========================================" << LogIO::POST;
609 0 : os << "No corrupting-type information has been set" << LogIO::POST;
610 0 : os << "===========================================" << LogIO::POST;
611 0 : return false;
612 : }
613 : else {
614 0 : os << "Visibilities will be CORRUPTED with the following terms:" << LogIO::POST;
615 :
616 0 : Int napp(vc_p.nelements());
617 0 : for (Int iapp=0;iapp<napp;++iapp)
618 : os << LogIO::NORMAL << ". "
619 0 : << vc_p[iapp]->siminfo()
620 0 : << LogIO::POST;
621 :
622 : // Report also on the noise settings
623 0 : noiseSummary(os);
624 :
625 : }
626 0 : return true;
627 : }
628 :
629 :
630 0 : Bool Simulator::noiseSummary(LogIO& os)
631 : {
632 0 : if (!ac_p) {
633 0 : return false;
634 : } else {
635 0 : os << "Thermal noise corruption activated" << LogIO::POST;
636 0 : os << "Thermal noise mode: " << noisemode_p << LogIO::POST;
637 : }
638 0 : return true;
639 : }
640 :
641 :
642 :
643 :
644 :
645 :
646 :
647 :
648 :
649 :
650 :
651 :
652 :
653 : //========================================================================
654 : // SETUP OBSERVATION
655 :
656 :
657 27 : Bool Simulator::settimes(const Quantity& integrationTime,
658 : const Bool useHourAngle,
659 : const MEpoch& refTime)
660 : {
661 :
662 81 : LogIO os(LogOrigin("simulator", "settimes()"));
663 : // Negative integration time crashes casa
664 27 : AlwaysAssert(integrationTime.getValue("s")>=0, AipsError);
665 : try {
666 :
667 27 : integrationTime_p=integrationTime;
668 27 : useHourAngle_p=useHourAngle;
669 27 : refTime_p=refTime;
670 :
671 : os << "Times " << endl
672 27 : << " Integration time " << integrationTime.getValue("s") << "s" << LogIO::POST;
673 27 : if(useHourAngle) {
674 27 : os << " Times will be interpreted as hour angles for first source" << LogIO::POST;
675 : }
676 :
677 27 : sim_p->settimes(integrationTime, useHourAngle, refTime);
678 :
679 27 : timesHaveBeenSet_p=true;
680 :
681 27 : return true;
682 0 : } catch (AipsError x) {
683 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
684 0 : return false;
685 : }
686 : return true;
687 :
688 : }
689 :
690 :
691 :
692 16 : Bool Simulator::setseed(const Int seed) {
693 16 : seed_p = seed;
694 16 : return true;
695 : }
696 :
697 :
698 :
699 28 : Bool Simulator::setconfig(const String& telname,
700 : const Vector<Double>& x,
701 : const Vector<Double>& y,
702 : const Vector<Double>& z,
703 : const Vector<Double>& dishDiameter,
704 : const Vector<Double>& offset,
705 : const Vector<String>& mount,
706 : const Vector<String>& antName,
707 : const Vector<String>& padName,
708 : const String& coordsystem,
709 : const MPosition& mRefLocation)
710 : {
711 56 : LogIO os(LogOrigin("Simulator", "setconfig()"));
712 :
713 28 : telescope_p = telname;
714 28 : x_p.resize(x.nelements());
715 28 : x_p = x;
716 28 : y_p.resize(y.nelements());
717 28 : y_p = y;
718 28 : z_p.resize(z.nelements());
719 28 : z_p = z;
720 28 : diam_p.resize(dishDiameter.nelements());
721 28 : diam_p = dishDiameter;
722 28 : offset_p.resize(offset.nelements());
723 28 : offset_p = offset;
724 28 : mount_p.resize(mount.nelements());
725 28 : mount_p = mount;
726 28 : antName_p.resize(antName.nelements());
727 28 : antName_p = antName;
728 28 : padName_p.resize(padName.nelements());
729 28 : padName_p = padName;
730 28 : coordsystem_p = coordsystem;
731 28 : mRefLocation_p = mRefLocation;
732 :
733 28 : uInt nn = x_p.nelements();
734 :
735 28 : if (diam_p.nelements() == 1) {
736 10 : diam_p.resize(nn);
737 10 : diam_p.set(dishDiameter(0));
738 : }
739 28 : if (mount_p.nelements() == 1) {
740 28 : mount_p.resize(nn);
741 28 : mount_p.set(mount(0));
742 : }
743 28 : if (mount_p.nelements() == 0) {
744 0 : mount_p.resize(nn);
745 0 : mount_p.set("alt-az");
746 : }
747 28 : if (offset_p.nelements() == 1) {
748 28 : offset_p.resize(nn);
749 28 : offset_p.set(offset(0));
750 : }
751 28 : if (offset_p.nelements() == 0) {
752 0 : offset_p.resize(nn);
753 0 : offset_p.set(0.0);
754 : }
755 28 : if (antName_p.nelements() == 1) {
756 10 : antName_p.resize(nn);
757 10 : antName_p.set(antName(0));
758 : }
759 28 : if (antName_p.nelements() == 0) {
760 0 : antName_p.resize(nn);
761 0 : antName_p.set("UNKNOWN");
762 : }
763 28 : if (padName_p.nelements() == 1) {
764 10 : padName_p.resize(nn);
765 10 : padName_p.set(padName(0));
766 : }
767 28 : if (padName_p.nelements() == 0) {
768 0 : padName_p.resize(nn);
769 0 : padName_p.set("UNKNOWN");
770 : }
771 :
772 28 : AlwaysAssert( (nn == y_p.nelements()) , AipsError);
773 28 : AlwaysAssert( (nn == z_p.nelements()) , AipsError);
774 28 : AlwaysAssert( (nn == diam_p.nelements()) , AipsError);
775 28 : AlwaysAssert( (nn == offset_p.nelements()) , AipsError);
776 28 : AlwaysAssert( (nn == mount_p.nelements()) , AipsError);
777 :
778 28 : areStationCoordsSet_p = true;
779 :
780 56 : sim_p->initAnt(telescope_p, x_p, y_p, z_p, diam_p, offset_p, mount_p, antName_p, padName_p,
781 28 : coordsystem_p, mRefLocation_p);
782 :
783 56 : return true;
784 : }
785 :
786 :
787 :
788 16 : Bool Simulator::getconfig() {
789 : // get it from NewMSSimulator
790 32 : Matrix<Double> xyz_p;
791 : Int nAnt;
792 16 : if (sim_p->getAnt(telescope_p, nAnt, &xyz_p, diam_p, offset_p, mount_p, antName_p, padName_p,
793 16 : coordsystem_p, mRefLocation_p)) {
794 16 : x_p.resize(nAnt);
795 16 : y_p.resize(nAnt);
796 16 : z_p.resize(nAnt);
797 285 : for (Int i=0;i<nAnt;i++) {
798 269 : x_p(i)=xyz_p(0,i);
799 269 : y_p(i)=xyz_p(1,i);
800 269 : z_p(i)=xyz_p(2,i);
801 : }
802 16 : areStationCoordsSet_p = true;
803 16 : return true;
804 : } else {
805 0 : return false;
806 : }
807 : }
808 :
809 :
810 :
811 488 : Bool Simulator::setfield(const String& sourceName,
812 : const MDirection& sourceDirection,
813 : const String& calCode,
814 : const Quantity& distance)
815 : {
816 1464 : LogIO os(LogOrigin("Simulator", "setfield()"));
817 : try {
818 :
819 488 : if (sourceName == "") {
820 0 : os << LogIO::SEVERE << "must provide a source name" << LogIO::POST;
821 0 : return false;
822 : }
823 :
824 488 : nField++;
825 :
826 488 : if (prtlev()>2) os << "nField = " << nField << LogIO::POST;
827 :
828 488 : distance_p.resize(nField,true);
829 488 : distance_p[nField-1]=distance;
830 488 : sourceName_p.resize(nField,true);
831 488 : sourceName_p[nField-1]=sourceName;
832 488 : sourceDirection_p.resize(nField,true);
833 488 : sourceDirection_p[nField-1]=sourceDirection;
834 488 : calCode_p.resize(nField,true);
835 488 : calCode_p[nField-1]=calCode;
836 :
837 488 : sim_p->initFields(sourceName, sourceDirection, calCode);
838 :
839 0 : } catch (AipsError x) {
840 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
841 0 : return false;
842 : }
843 488 : return true;
844 : };
845 :
846 :
847 :
848 :
849 0 : Bool Simulator::setmosaicfield(const String& sourcename, const String& calcode, const MDirection& fieldcenter, const Int xmosp, const Int ymosp, const Quantity& mosspacing, const Quantity& distance)
850 : {
851 :
852 0 : Int nx = xmosp;
853 0 : Int ny = ymosp;
854 0 : Double roffset = mosspacing.getValue("rad");
855 : Double newraval, newdecval;
856 0 : uInt k=0;
857 0 : for (Int i=0; i< nx; ++i) {
858 0 : for(Int j=0; j < ny; ++j) {
859 0 : if((nx/2)!=floor(nx/2)) { // odd number of fields in x direction(ra)
860 0 : newraval = fieldcenter.getAngle().getValue("rad")(0) +
861 0 : (i-ceil(nx/2.0))*roffset/cos(fieldcenter.getAngle().getValue("rad")(1));
862 : }
863 : else { //even case
864 0 : newraval = fieldcenter.getAngle().getValue("rad")(0) +
865 0 : ((i-ceil(nx/2)) - 0.5)*roffset/cos(fieldcenter.getAngle().getValue("rad")(1));
866 : }
867 0 : if((ny/2)!=floor(ny/2)) {
868 0 : newdecval = fieldcenter.getAngle().getValue("rad")(1) +
869 0 : (j-ceil(ny/2.0))*roffset;
870 : }
871 : else {
872 0 : newdecval = fieldcenter.getAngle().getValue("rad")(1) +
873 0 : ((j-ceil(ny/2.0)) - 0.5)*roffset;
874 : }
875 0 : if(newraval >2.0*C::pi) {
876 0 : newraval = newraval - 2.0*C::pi;
877 : }
878 : Int sign;
879 0 : if(abs(newdecval) >C::pi/2) {
880 0 : if(newdecval<0) {
881 0 : sign = -1;
882 : }
883 : else {
884 0 : sign = 1;
885 : }
886 0 : newdecval = sign*(C::pi - abs(newdecval));
887 0 : newraval = abs(C::pi - newraval);
888 : }
889 0 : Quantity newdirra(newraval, "rad");
890 0 : Quantity newdirdec(newdecval, "rad");
891 0 : MDirection newdir(newdirra, newdirdec);
892 0 : newdir.setRefString(fieldcenter.getRefString());
893 0 : ostringstream oos;
894 0 : oos << sourcename << "_" << k ;
895 :
896 :
897 0 : setfield(String(oos), newdir, calcode, distance);
898 :
899 0 : ++k;
900 : }
901 : }
902 :
903 :
904 0 : return true;
905 : }
906 :
907 :
908 :
909 28 : Bool Simulator::setspwindow(const String& spwName,
910 : const Quantity& freq,
911 : const Quantity& deltafreq,
912 : const Quantity& freqresolution,
913 : const MFrequency::Types& freqType,
914 : const Int nChan,
915 : const String& stokes)
916 :
917 : {
918 84 : LogIO os(LogOrigin("Simulator", "setspwindow()"));
919 : try {
920 28 : if (nChan == 0) {
921 0 : os << LogIO::SEVERE << "must provide nchannels" << LogIO::POST;
922 0 : return false;
923 : }
924 :
925 28 : nSpw++;
926 : #ifdef RI_DEBUG
927 : os << "nspw = " << nSpw << LogIO::POST;
928 : #endif
929 28 : spWindowName_p.resize(nSpw,true);
930 28 : spWindowName_p[nSpw-1] = spwName;
931 28 : nChan_p.resize(nSpw,true);
932 28 : nChan_p[nSpw-1] = nChan;
933 28 : startFreq_p.resize(nSpw,true);
934 28 : startFreq_p[nSpw-1] = freq;
935 28 : freqInc_p.resize(nSpw,true);
936 28 : freqInc_p[nSpw-1] = deltafreq;
937 28 : freqRes_p.resize(nSpw,true);
938 28 : freqRes_p[nSpw-1] = freqresolution;
939 28 : stokesString_p.resize(nSpw,true);
940 28 : stokesString_p[nSpw-1] = stokes;
941 :
942 : #ifdef RI_DEBUG
943 : os << "sending init to MSSim for spw = " << spWindowName_p[nSpw-1] << LogIO::POST;
944 : #endif
945 :
946 : //freqType=MFrequency::TOPO;
947 28 : sim_p->initSpWindows(spWindowName_p[nSpw-1], nChan_p[nSpw-1],
948 28 : startFreq_p[nSpw-1], freqInc_p[nSpw-1],
949 28 : freqRes_p[nSpw-1], freqType,
950 28 : stokesString_p[nSpw-1]);
951 :
952 0 : } catch (AipsError x) {
953 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
954 0 : return false;
955 : }
956 28 : return true;
957 : };
958 :
959 :
960 28 : Bool Simulator::setfeed(const String& mode,
961 : const Vector<Double>& x,
962 : const Vector<Double>& y,
963 : const Vector<String>& pol)
964 : {
965 84 : LogIO os(LogOrigin("Simulator", "setfeed()"));
966 :
967 28 : if (mode != "perfect R L" && mode != "perfect X Y" && mode != "list") {
968 : os << LogIO::SEVERE <<
969 : "Currently, only perfect R L or perfect X Y feeds or list are recognized"
970 0 : << LogIO::POST;
971 0 : return false;
972 : }
973 28 : feedMode_p = mode;
974 28 : sim_p->initFeeds(feedMode_p, x, y, pol);
975 :
976 28 : nFeeds_p = x.nelements();
977 28 : feedsHaveBeenSet = true;
978 :
979 28 : return true;
980 : };
981 :
982 :
983 :
984 :
985 17 : Bool Simulator::setvp(const Bool dovp,
986 : const Bool doDefaultVPs,
987 : const String& vpTable,
988 : const Bool doSquint,
989 : const Quantity &parAngleInc,
990 : const Quantity &skyPosThreshold,
991 : const Float &pbLimit)
992 : {
993 34 : LogIO os(LogOrigin("Simulator", "setvp()"));
994 :
995 17 : os << "Setting voltage pattern parameters" << LogIO::POST;
996 17 : doVP_p=dovp;
997 17 : doDefaultVP_p = doDefaultVPs;
998 17 : vpTableStr_p = vpTable;
999 17 : if (doSquint) {
1000 17 : squintType_p = BeamSquint::GOFIGURE;
1001 : } else {
1002 0 : squintType_p = BeamSquint::NONE;
1003 : }
1004 17 : parAngleInc_p = parAngleInc;
1005 17 : skyPosThreshold_p = skyPosThreshold;
1006 :
1007 17 : if (doDefaultVP_p) {
1008 0 : os << "Using system default voltage patterns for each telescope" << LogIO::POST;
1009 : } else {
1010 17 : if (vpTableStr_p != String("")) {
1011 0 : os << "Using user defined voltage patterns in Table "<< vpTableStr_p
1012 0 : << LogIO::POST;
1013 : }
1014 : }
1015 17 : if (doSquint) {
1016 17 : os << "Beam Squint will be included in the VP model" << LogIO::POST;
1017 : os << "and the parallactic angle increment is "
1018 17 : << parAngleInc_p.getValue("deg") << " degrees" << LogIO::POST;
1019 : }
1020 17 : pbLimit_p = pbLimit;
1021 34 : return true;
1022 : };
1023 :
1024 :
1025 :
1026 :
1027 :
1028 :
1029 :
1030 :
1031 :
1032 :
1033 :
1034 :
1035 :
1036 : //========================================================================
1037 : // Create corruption terms - top level specialized routines
1038 :
1039 :
1040 : // NEW NOISE WITH ANoise
1041 :
1042 16 : Bool Simulator::setnoise(const String& mode,
1043 : const String& caltable, // if blank, not stored
1044 : const Quantity& simplenoise,
1045 : // ATM calculation
1046 : const Quantity& pground,
1047 : const Float relhum,
1048 : const Quantity& altitude,
1049 : const Quantity& waterheight,
1050 : const Quantity& pwv,
1051 : // OR user-specified tau and tatmos
1052 : const Float tatmos=250.0,
1053 : const Float tau=0.1,
1054 : // antenna parameters used for either option
1055 : const Float antefficiency=0.80,
1056 : const Float spillefficiency=0.85,
1057 : const Float correfficiency=0.85,
1058 : const Float trx=150.0,
1059 : const Float tground=270.0,
1060 : const Float tcmb=2.73,
1061 : const Bool OTF=true,
1062 : const Float senscoeff=0.0,
1063 : const Int rxtype=0
1064 : ) {
1065 :
1066 48 : LogIO os(LogOrigin("Simulator", "setnoise2()", WHERE));
1067 :
1068 : try {
1069 :
1070 : //cout<<" Sim::setnoise "<<pground<<" "<<relhum<<" "<<altitude<<" "<<waterheight<<" "<<pwv<<endl;
1071 :
1072 16 : noisemode_p = mode;
1073 :
1074 32 : RecordDesc simparDesc;
1075 16 : simparDesc.addField ("type", TpString);
1076 16 : simparDesc.addField ("mode", TpString);
1077 16 : simparDesc.addField ("caltable", TpString);
1078 :
1079 16 : simparDesc.addField ("amplitude", TpFloat); // for constant scale
1080 : // simparDesc.addField ("scale", TpFloat); // for fractional fluctuations
1081 16 : simparDesc.addField ("combine" ,TpString);
1082 :
1083 : // have to be defined here or else have to be added later
1084 16 : simparDesc.addField ("startTime", TpDouble);
1085 16 : simparDesc.addField ("stopTime", TpDouble);
1086 :
1087 16 : simparDesc.addField ("antefficiency" ,TpFloat);
1088 16 : simparDesc.addField ("spillefficiency",TpFloat);
1089 16 : simparDesc.addField ("correfficiency" ,TpFloat);
1090 16 : simparDesc.addField ("trx" ,TpFloat);
1091 16 : simparDesc.addField ("tground" ,TpFloat);
1092 16 : simparDesc.addField ("tcmb" ,TpFloat);
1093 16 : simparDesc.addField ("senscoeff" ,TpFloat);
1094 16 : simparDesc.addField ("rxType" ,TpInt);
1095 :
1096 : // user-override of ATM calculated tau
1097 16 : simparDesc.addField ("tatmos" ,TpFloat);
1098 16 : simparDesc.addField ("tau0" ,TpFloat);
1099 :
1100 16 : simparDesc.addField ("mean_pwv" ,TpDouble);
1101 16 : simparDesc.addField ("pground" ,TpDouble);
1102 16 : simparDesc.addField ("relhum" ,TpFloat);
1103 16 : simparDesc.addField ("altitude" ,TpDouble);
1104 16 : simparDesc.addField ("waterheight" ,TpDouble);
1105 :
1106 16 : simparDesc.addField ("seed" ,TpInt);
1107 :
1108 : // RI todo setnoise2 if tau0 is not defined, use freqdep
1109 :
1110 32 : String caltbl=caltable;
1111 16 : caltbl.trim();
1112 : string::size_type strlen;
1113 16 : strlen=caltbl.length();
1114 16 : if (strlen>3)
1115 10 : if (caltbl.substr(strlen-3,3)=="cal") {
1116 0 : caltbl.resize(strlen-3);
1117 0 : strlen-=3;
1118 : }
1119 16 : if (strlen>1)
1120 10 : if (caltbl.substr(strlen-1,1)==".") {
1121 0 : caltbl.resize(strlen-1);
1122 0 : strlen-=1;
1123 : }
1124 16 : if (strlen>1)
1125 10 : if (caltbl.substr(strlen-1,1)=="_") {
1126 0 : caltbl.resize(strlen-1);
1127 0 : strlen-=1;
1128 : }
1129 :
1130 32 : Record simpar(simparDesc,RecordInterface::Variable);
1131 16 : simpar.define ("seed", seed_p);
1132 16 : simpar.define ("type", "A Noise");
1133 16 : if (strlen>1)
1134 10 : simpar.define ("caltable", caltbl+".A.cal");
1135 16 : simpar.define ("mode", mode);
1136 16 : simpar.define ("combine", ""); // SPW,FIELD, etc
1137 :
1138 16 : if (mode=="simplenoise") {
1139 : os << "Using simple noise model with noise level of " << simplenoise.getValue("Jy")
1140 0 : << " Jy" << LogIO::POST;
1141 0 : simpar.define ("amplitude", Float(simplenoise.getValue("Jy")) );
1142 0 : simpar.define ("mode", "simple");
1143 :
1144 16 : } else if (mode=="tsys-atm" or mode=="tsys-manual") {
1145 : // do be scaled in a minute by a Tsys-derived M below
1146 16 : simpar.define ("amplitude", Float(1.0) );
1147 : // os << "adding noise with unity amplitude" << LogIO::POST;
1148 16 : simpar.define ("mode", mode);
1149 :
1150 : } else {
1151 0 : throw(AipsError("unsupported mode "+mode+" in setnoise()"));
1152 : }
1153 :
1154 16 : Bool saveOnthefly=false;
1155 16 : if (simpar.isDefined("onthefly")) {
1156 0 : saveOnthefly=simpar.asBool("onthefly");
1157 : }
1158 16 : simpar.define("onthefly",OTF);
1159 :
1160 : // create the ANoise
1161 16 : if (!create_corrupt(simpar))
1162 0 : throw(AipsError("could not create ANoise in Simulator::setnoise"));
1163 :
1164 16 : if (mode=="tsys-atm" or mode=="tsys-manual") {
1165 :
1166 16 : simpar.define ("antefficiency" ,antefficiency );
1167 16 : simpar.define ("correfficiency" ,correfficiency );
1168 16 : simpar.define ("spillefficiency",spillefficiency);
1169 16 : simpar.define ("trx" ,trx );
1170 16 : simpar.define ("tground" ,tground );
1171 16 : simpar.define ("tcmb" ,tcmb );
1172 :
1173 16 : if (rxtype>=0) {
1174 16 : simpar.define ("rxType", rxtype);
1175 16 : if (rxtype>0) {
1176 0 : os<<"User has requested Double Sideband Receiver"<<LogIO::POST;
1177 : }
1178 : } else {
1179 0 : simpar.define ("rxType", 0);
1180 0 : os<<"User has not set Rx type, using 2SB"<<LogIO::POST;
1181 : }
1182 :
1183 16 : if ( senscoeff > 0.0 ) {
1184 10 : simpar.define ("senscoeff", Float(senscoeff) );
1185 10 : os << "adding noise with the sensitivity constant of " << senscoeff << LogIO::POST;
1186 : } else {
1187 6 : simpar.define("senscoeff", Float(1./ sqrt(2.)));
1188 6 : os << "adding noise with the sensitivity constant of 1/sqrt(2)" << LogIO::POST;
1189 : }
1190 :
1191 16 : if (pwv.getValue("mm")>0.) {
1192 10 : if (pwv.getValue("mm")>100.)
1193 0 : throw(AipsError(" you input PWV="+String::toString(pwv.getValue("mm"))+"mm. The most common reason for this error is forgetting to set units, which default to meters"));
1194 :
1195 10 : simpar.define ("mean_pwv", pwv.getValue("mm"));
1196 : } else {
1197 6 : simpar.define ("mean_pwv", Double(1.));
1198 : // we want to set it, but it doesn't get used unless ATM is being used
1199 6 : if (mode=="tsys-atm") os<<"User has not set PWV, using 1mm"<<LogIO::POST;
1200 : }
1201 :
1202 16 : if (mode=="tsys-manual") {
1203 : // user can override the ATM calculated optical depth
1204 : // with tau0 to be used over the entire SPW,
1205 6 : simpar.define ("tau0" ,tau );
1206 6 : if (tatmos>10)
1207 6 : simpar.define ("tatmos" ,tatmos );
1208 : else
1209 0 : simpar.define ("tatmos" ,250. );
1210 : // AtmosCorruptor cal deal with
1211 : // an MF in tsys-manual mode - it will use ATM to calculate
1212 : // the relative opacity across the band, but it won't properly
1213 : // integrate the atmosphere to get T_ebb.
1214 : // so for now we'll just make tsys-manual mean freqDepPar=false
1215 :
1216 6 : simpar.define ("type", "T");
1217 : //simpar.define ("type", "M");
1218 :
1219 : } else {
1220 : // otherwise ATM will be used to calculate tau from pwv
1221 : // catch uninitialized variant @#$^@! XML interface
1222 10 : if (pground.getValue("mbar")>100.)
1223 0 : simpar.define ("pground", pground.getValue("mbar"));
1224 : else {
1225 10 : simpar.define ("pground", Double(560.));
1226 10 : os<<"User has not set ground pressure, using 560mb"<<LogIO::POST;
1227 : }
1228 10 : if (relhum>0)
1229 10 : simpar.define ("relhum", relhum);
1230 : else {
1231 0 : simpar.define ("relhum", 20.);
1232 0 : os<<"User has not set ground relative humidity, using 20%"<<LogIO::POST;
1233 : }
1234 10 : if (altitude.getValue("m")>0.)
1235 0 : simpar.define ("altitude", altitude.getValue("m"));
1236 : else {
1237 10 : simpar.define ("altitude", Double(5000.));
1238 10 : os<<"User has not set site altitude, using 5000m"<<LogIO::POST;
1239 : }
1240 10 : if (waterheight.getValue("m")>100.)
1241 0 : simpar.define ("waterheight", waterheight.getValue("km"));
1242 : else {
1243 10 : simpar.define ("waterheight", Double(2.0));
1244 10 : os<<"User has not set water scale height, using 2km"<<LogIO::POST;
1245 : }
1246 : // as a function of frequency (freqDepPar=true)
1247 : //simpar.define ("type", "TF");
1248 10 : simpar.define ("type", "TF NOISE");
1249 : // simpar.define ("type", "MF");
1250 : }
1251 :
1252 16 : if (strlen>1)
1253 10 : simpar.define ("caltable", caltbl+".T.cal");
1254 : //simpar.define ("caltable", caltbl+".M.cal");
1255 :
1256 16 : simpar.define("onthefly",saveOnthefly);
1257 :
1258 : // create the T
1259 16 : if (!create_corrupt(simpar))
1260 0 : throw(AipsError("could not create T in Simulator::setnoise"));
1261 : //throw(AipsError("could not create M in Simulator::setnoise"));
1262 : }
1263 :
1264 0 : } catch (AipsError x) {
1265 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1266 0 : return false;
1267 : }
1268 16 : return true;
1269 : }
1270 :
1271 :
1272 :
1273 :
1274 0 : Bool Simulator::setgain(const String& mode,
1275 : const String& caltable,
1276 : const Quantity& interval,
1277 : const Vector<Double>& amplitude)
1278 : {
1279 0 : LogIO os(LogOrigin("Simulator", "setgain()", WHERE));
1280 :
1281 : try {
1282 :
1283 0 : if(mode=="table") {
1284 0 : os << LogIO::SEVERE << "Cannot yet read from table" << LogIO::POST;
1285 0 : return false;
1286 : }
1287 : else {
1288 : // RI TODO Sim::setgain add mode=simple and =normal
1289 0 : if (mode=="fbm" or mode=="random") {
1290 :
1291 : // set record format for calibration table simulation information
1292 0 : RecordDesc simparDesc;
1293 0 : simparDesc.addField ("type", TpString);
1294 0 : simparDesc.addField ("caltable", TpString);
1295 0 : simparDesc.addField ("mode", TpString);
1296 0 : simparDesc.addField ("interval", TpDouble);
1297 0 : simparDesc.addField ("camp", TpComplex);
1298 0 : simparDesc.addField ("amplitude", TpDouble);
1299 0 : simparDesc.addField ("combine", TpString);
1300 0 : simparDesc.addField ("startTime", TpDouble);
1301 0 : simparDesc.addField ("stopTime", TpDouble);
1302 0 : simparDesc.addField ("seed", TpInt);
1303 :
1304 : // Create record with the requisite field values
1305 0 : Record simpar(simparDesc,RecordInterface::Variable);
1306 0 : simpar.define ("type", "G JONES");
1307 0 : simpar.define ("interval", interval.getValue("s"));
1308 0 : simpar.define ("mode", mode);
1309 0 : Complex camp(0.1,0.1);
1310 0 : if (amplitude.size()==1)
1311 0 : camp=Complex(amplitude[0],amplitude[0]);
1312 : else
1313 0 : camp=Complex(amplitude[0],amplitude[1]);
1314 0 : simpar.define ("camp", camp);
1315 0 : os << LogIO::NORMAL << "Gain corruption with complex RMS amplitude = " << camp << LogIO::POST;
1316 0 : simpar.define ("amplitude", Double(abs(camp)));
1317 : //simpar.define ("amplitude", amplitude);
1318 0 : simpar.define ("caltable", caltable);
1319 0 : simpar.define ("combine", "");
1320 0 : simpar.define ("seed", seed_p);
1321 :
1322 : // create the G
1323 0 : if (!create_corrupt(simpar))
1324 0 : throw(AipsError("could not create G in Simulator::setgain"));
1325 :
1326 : } else {
1327 0 : throw(AipsError("unsupported mode "+mode+" in setgain()"));
1328 : }
1329 : }
1330 0 : } catch (AipsError x) {
1331 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1332 0 : return false;
1333 : }
1334 0 : return true;
1335 : }
1336 :
1337 :
1338 :
1339 :
1340 :
1341 :
1342 :
1343 0 : Bool Simulator::settrop(const String& mode,
1344 : const String& caltable, // output
1345 : const Float pwv,
1346 : const Float deltapwv,
1347 : const Float beta,
1348 : const Float windspeed,
1349 : const Float simintsec=-1.) {
1350 :
1351 0 : LogIO os(LogOrigin("Simulator", "settrop()", WHERE));
1352 :
1353 : #ifndef RI_DEBUG
1354 : try {
1355 : #endif
1356 :
1357 0 : if (mode=="test"||mode=="individual"||mode=="screen") {
1358 :
1359 : // set record format for calibration table simulation information
1360 0 : RecordDesc simparDesc;
1361 0 : simparDesc.addField ("type", TpString);
1362 0 : simparDesc.addField ("caltable", TpString);
1363 0 : simparDesc.addField ("mean_pwv", TpFloat);
1364 0 : simparDesc.addField ("mode", TpString);
1365 0 : simparDesc.addField ("delta_pwv", TpFloat);
1366 0 : simparDesc.addField ("beta", TpFloat);
1367 0 : simparDesc.addField ("windspeed", TpFloat);
1368 0 : simparDesc.addField ("combine", TpString);
1369 :
1370 0 : if(simintsec>0.){
1371 0 : simparDesc.addField ("simint", TpString);
1372 : }
1373 :
1374 0 : simparDesc.addField ("startTime", TpDouble);
1375 0 : simparDesc.addField ("stopTime", TpDouble);
1376 :
1377 0 : simparDesc.addField ("antefficiency" ,TpFloat);
1378 0 : simparDesc.addField ("spillefficiency",TpFloat);
1379 0 : simparDesc.addField ("correfficiency" ,TpFloat);
1380 0 : simparDesc.addField ("trx" ,TpFloat);
1381 0 : simparDesc.addField ("tcmb" ,TpFloat);
1382 0 : simparDesc.addField ("tatmos" ,TpFloat);
1383 :
1384 0 : simparDesc.addField ("tground" ,TpFloat);
1385 0 : simparDesc.addField ("pground" ,TpDouble);
1386 0 : simparDesc.addField ("relhum" ,TpFloat);
1387 0 : simparDesc.addField ("altitude" ,TpDouble);
1388 0 : simparDesc.addField ("waterheight" ,TpDouble);
1389 :
1390 0 : simparDesc.addField ("seed" ,TpInt);
1391 :
1392 : // create record with the requisite field values
1393 0 : Record simpar(simparDesc,RecordInterface::Variable);
1394 0 : simpar.define ("type", "TF");
1395 0 : simpar.define ("caltable", caltable);
1396 0 : simpar.define ("mean_pwv", pwv);
1397 0 : simpar.define ("mode", mode);
1398 0 : simpar.define ("delta_pwv", deltapwv);
1399 0 : simpar.define ("beta", beta);
1400 0 : simpar.define ("windspeed", windspeed);
1401 0 : simpar.define ("combine", "");
1402 :
1403 0 : if(simintsec>0.){
1404 0 : simpar.define ("simint", String::toString(simintsec)+"s");
1405 : }
1406 :
1407 0 : simpar.define ("seed", seed_p);
1408 :
1409 : // if (tground>100.)
1410 : // simpar.define ("tground", tground);
1411 : // else {
1412 0 : simpar.define ("tground", Float(270.));
1413 : // os<<"User has not set ground temperature, using 270K"<<LogIO::POST;
1414 : // }
1415 : // if (pground.getValue("mbar")>100.)
1416 : // simpar.define ("pground", pground.getValue("mbar"));
1417 : // else {
1418 0 : simpar.define ("pground", Double(560.));
1419 : // os<<"User has not set ground pressure, using 560mb"<<LogIO::POST;
1420 : // }
1421 : // if (relhum>0)
1422 : // simpar.define ("relhum", relhum);
1423 : // else {
1424 0 : simpar.define ("relhum", Float(20.));
1425 : // os<<"User has not set ground relative humidity, using 20%"<<LogIO::POST;
1426 : // }
1427 : // if (altitude.getValue("m")>0.)
1428 : // simpar.define ("altitude", altitude.getValue("m"));
1429 : // else {
1430 0 : simpar.define ("altitude", Double(5000.));
1431 : // os<<"User has not set site altitude, using 5000m"<<LogIO::POST;
1432 : // }
1433 : // if (waterheight.getValue("m")>100.)
1434 : // simpar.define ("waterheight", waterheight.getValue("km"));
1435 : // else {
1436 0 : simpar.define ("waterheight", Double(2.0)); // km
1437 : // os<<"User has not set water scale height, using 2km"<<LogIO::POST;
1438 : // }
1439 0 : simpar.define ("spillefficiency", Float(.85));
1440 :
1441 : // create the T
1442 0 : if (!create_corrupt(simpar))
1443 0 : throw(AipsError("could not create T in Simulator::settrop"));
1444 :
1445 : } else {
1446 0 : throw(AipsError("unsupported mode "+mode+" in settrop()"));
1447 : }
1448 :
1449 : #ifndef RI_DEBUG
1450 0 : } catch (AipsError x) {
1451 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1452 0 : return false;
1453 : }
1454 : #endif
1455 0 : return true;
1456 : }
1457 :
1458 :
1459 :
1460 :
1461 :
1462 :
1463 0 : Bool Simulator::setleakage(const String& /*mode*/, const String& table,
1464 : //const Quantity& interval,
1465 : const Vector<Double>& amplitude,
1466 : const Vector<Double>& offset)
1467 : {
1468 :
1469 0 : LogIO os(LogOrigin("Simulator", "setleakage()", WHERE));
1470 :
1471 : #ifndef RI_DEBUG
1472 : try {
1473 : #endif
1474 :
1475 : // set record format for calibration table simulation information
1476 0 : RecordDesc simparDesc;
1477 0 : simparDesc.addField ("type", TpString);
1478 0 : simparDesc.addField ("caltable", TpString);
1479 : // leave this one for generic SVC::createCorruptor
1480 0 : simparDesc.addField ("amplitude", TpDouble);
1481 0 : simparDesc.addField ("camp", TpComplex);
1482 0 : simparDesc.addField ("offset", TpComplex);
1483 0 : simparDesc.addField ("combine", TpString);
1484 : //simparDesc.addField ("interval", TpDouble);
1485 0 : simparDesc.addField ("simint", TpString);
1486 0 : simparDesc.addField ("startTime", TpDouble);
1487 0 : simparDesc.addField ("stopTime", TpDouble);
1488 :
1489 0 : simparDesc.addField ("seed", TpInt);
1490 :
1491 : // create record with the requisite field values
1492 0 : Record simpar(simparDesc,RecordInterface::Variable);
1493 0 : simpar.define ("type", "D");
1494 0 : simpar.define ("caltable", table);
1495 0 : Complex camp(0.1,0.1);
1496 0 : if (amplitude.size()==1)
1497 0 : camp=Complex(amplitude[0],amplitude[0]);
1498 : else
1499 0 : camp=Complex(amplitude[0],amplitude[1]);
1500 0 : simpar.define ("camp", camp);
1501 0 : simpar.define ("amplitude", Double(abs(camp)));
1502 0 : Complex off(0.,0.);
1503 0 : if (offset.size()==1)
1504 0 : off=Complex(offset[0],offset[0]);
1505 : else
1506 0 : off=Complex(offset[0],offset[1]);
1507 0 : simpar.define ("offset", off);
1508 :
1509 : //simpar.define ("interval", interval.getValue("s"));
1510 : // provide infinite interval
1511 0 : simpar.define ("simint", "infinite");
1512 :
1513 0 : simpar.define ("combine", "");
1514 0 : simpar.define ("seed", seed_p);
1515 :
1516 :
1517 : // create the D
1518 0 : if (!create_corrupt(simpar))
1519 0 : throw(AipsError("could not create D in Simulator::setleakage"));
1520 :
1521 : #ifndef RI_DEBUG
1522 0 : } catch (AipsError x) {
1523 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1524 0 : return false;
1525 : }
1526 : #endif
1527 0 : return true;
1528 : }
1529 :
1530 :
1531 :
1532 :
1533 :
1534 :
1535 :
1536 :
1537 :
1538 :
1539 :
1540 :
1541 :
1542 :
1543 :
1544 :
1545 :
1546 :
1547 :
1548 : //========================================================================
1549 : // OLD as yet unconverted corruption generation routines
1550 :
1551 : // OLD NOISE WITH ACoh
1552 :
1553 0 : Bool Simulator::oldsetnoise(const String& mode,
1554 : const String& /*table*/,
1555 : const Quantity& simplenoise,
1556 : const Float antefficiency=0.80,
1557 : const Float correfficiency=0.85,
1558 : const Float spillefficiency=0.85,
1559 : const Float tau=0.0,
1560 : const Float trx=50.0,
1561 : const Float tatmos=250.0,
1562 : const Float tcmb=2.73) {
1563 :
1564 0 : LogIO os(LogOrigin("Simulator", "oldsetnoise()", WHERE));
1565 : try {
1566 :
1567 0 : noisemode_p = mode;
1568 :
1569 0 : os << LogIO::WARN << "Using deprecated ACoh Noise - this will dissapear in the future - please switch to sm.setnoise unless you are simulating single dish data" << LogIO::POST;
1570 :
1571 0 : if(mode=="table") {
1572 0 : os << LogIO::SEVERE << "Cannot yet read from table" << LogIO::POST;
1573 0 : return false;
1574 : }
1575 0 : else if (mode=="simplenoise") {
1576 : os << "Using simple noise model with noise level of " << simplenoise.getValue("Jy")
1577 0 : << " Jy" << LogIO::POST;
1578 0 : if(ac_p) delete ac_p; ac_p = 0;
1579 0 : ac_p = new SimACoh(seed_p, simplenoise.getValue("Jy") );
1580 : }
1581 : else {
1582 0 : os << "Using the Brown calculated noise model" << LogIO::POST;
1583 0 : os << " eta_ant=" << antefficiency << " eta_corr=" << correfficiency << " eta_spill=" << spillefficiency << LogIO::POST;
1584 0 : os << " tau=" << tau << " trx=" << trx << " tatmos=" << tatmos << " tcmb=" << tcmb << LogIO::POST;
1585 0 : if(ac_p) delete ac_p; ac_p = 0;
1586 0 : ac_p = new SimACohCalc(seed_p, antefficiency, correfficiency,
1587 0 : spillefficiency, tau, Quantity(trx, "K"),
1588 0 : Quantity(tatmos, "K"), Quantity(tcmb, "K"));
1589 : }
1590 :
1591 0 : return true;
1592 0 : } catch (AipsError x) {
1593 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1594 0 : return false;
1595 : }
1596 : return true;
1597 :
1598 : }
1599 :
1600 :
1601 :
1602 :
1603 :
1604 0 : Bool Simulator::setpa(const String& /*mode*/, const String& /*table*/,
1605 : const Quantity& /*interval*/) {
1606 :
1607 0 : LogIO os(LogOrigin("Simulator", "setpa()", WHERE));
1608 :
1609 : try {
1610 :
1611 0 : throw(AipsError("Corruption by simulated errors temporarily disabled (06Nov20 gmoellen)"));
1612 : /*
1613 : if(mode=="table") {
1614 : os << LogIO::SEVERE << "Cannot yet read from table" << LogIO::POST;
1615 : return false;
1616 : }
1617 : else {
1618 : makeVisSet();
1619 : if(pj_p) delete pj_p; pj_p = 0;
1620 : pj_p = new PJones (*vs_p, interval.get("s").getValue());
1621 : os <<"Using parallactic angle correction"<< LogIO::POST;
1622 : }
1623 : */
1624 : return true;
1625 0 : } catch (AipsError x) {
1626 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1627 0 : return false;
1628 : }
1629 : return true;
1630 : };
1631 :
1632 :
1633 :
1634 :
1635 0 : Bool Simulator::setbandpass(const String& /*mode*/, const String& /*table*/,
1636 : const Quantity& /*interval*/,
1637 : const Vector<Double>& /*amplitude*/) {
1638 :
1639 0 : LogIO os(LogOrigin("Simulator", "setbandpass()", WHERE));
1640 :
1641 : try {
1642 :
1643 0 : throw(AipsError("Corruption by simulated errors temporarily disabled (06Nov20 gmoellen)"));
1644 :
1645 : /*
1646 : if(mode=="table") {
1647 : os << LogIO::SEVERE << "Cannot yet read from table" << LogIO::POST;
1648 : return false;
1649 : }
1650 : else {
1651 : os << LogIO::SEVERE << "Cannot yet calculate bandpass" << LogIO::POST;
1652 : return false;
1653 : }
1654 : return true;
1655 : */
1656 0 : } catch (AipsError x) {
1657 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1658 :
1659 0 : return false;
1660 : }
1661 : return true;
1662 : }
1663 :
1664 :
1665 :
1666 :
1667 0 : Bool Simulator::setpointingerror(const String& epJTableName,
1668 : const Bool applyPointingOffsets,
1669 : const Bool doPBCorrection)
1670 : {
1671 0 : LogIO os(LogOrigin("Simulator", "close()", WHERE));
1672 0 : epJTableName_p = epJTableName;
1673 : // makeVisSet();
1674 : try {
1675 0 : throw(AipsError("Corruption by simulated errors temporarily disabled (06Nov20 gmoellen)"));
1676 : /*
1677 : if (epJ_p) delete epJ_p;epJ_p=0;
1678 : {
1679 : epJ_p = new EPJones(*vs_p);
1680 : epJ_p->load(epJTableName_p,"","diagonal");
1681 : }
1682 : */
1683 : }
1684 0 : catch (AipsError x)
1685 : {
1686 : os << LogIO::SEVERE << "Caught exception: "
1687 0 : << x.getMesg() << LogIO::POST;
1688 0 : return false;
1689 : }
1690 :
1691 : applyPointingOffsets_p = applyPointingOffsets;
1692 : doPBCorrection_p = doPBCorrection;
1693 : return true;
1694 : }
1695 :
1696 :
1697 :
1698 :
1699 :
1700 :
1701 :
1702 :
1703 :
1704 :
1705 :
1706 : //========================================================================
1707 : // CORRUPTION - GENERIC VISCAL FUNCTIONS
1708 :
1709 :
1710 :
1711 32 : Bool Simulator::create_corrupt(Record& simpar)
1712 : {
1713 64 : LogIO os(LogOrigin("Simulator", "create_corrupt()", WHERE));
1714 32 : SolvableVisCal *svc(NULL);
1715 :
1716 : // RI todo sim::create_corrupt assert that ms has certain structure
1717 :
1718 : try {
1719 32 : makeVisSet();
1720 :
1721 64 : String upType=simpar.asString("type");
1722 32 : upType.upcase();
1723 32 : os << LogIO::NORMAL << "Creating "<< upType <<" Calibration structure for data corruption." << LogIO::POST;
1724 :
1725 32 : svc = createSolvableVisCal(upType,*vs_p);
1726 :
1727 32 : svc->setPrtlev(prtlev());
1728 :
1729 64 : Vector<Double> solTimes;
1730 32 : svc->setSimulate(*vs_p,simpar,solTimes);
1731 :
1732 : // add to the pointer block of VCs:
1733 32 : uInt napp=vc_p.nelements();
1734 32 : vc_p.resize(napp+1,false,true);
1735 32 : vc_p[napp] = (VisCal*) svc;
1736 : // svc=NULL;
1737 32 : ve_p.setapply(vc_p);
1738 :
1739 0 : } catch (AipsError x) {
1740 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1741 0 : if (svc) delete svc;
1742 0 : throw(AipsError("Error in Simulator::create_corrupt"));
1743 : return false;
1744 : }
1745 32 : return true;
1746 : }
1747 :
1748 :
1749 :
1750 :
1751 :
1752 :
1753 :
1754 : //========================================================================
1755 : // corrupt and setapply, for actually changing visibilities
1756 :
1757 :
1758 : /// this can be used to load any table, just that it has to have the right form
1759 :
1760 0 : Bool Simulator::setapply(const String& type,
1761 : const Double& t,
1762 : const String& table,
1763 : const String& /*spw*/,
1764 : const String& /*field*/,
1765 : const String& interp,
1766 : const Bool& calwt,
1767 : const Vector<Int>& spwmap,
1768 : const Float& opacity)
1769 : {
1770 0 : LogIO os(LogOrigin("Simulator", "setapply()", WHERE));
1771 :
1772 : // First try to create the requested VisCal object
1773 0 : VisCal *vc(NULL);
1774 :
1775 : try {
1776 :
1777 0 : makeVisSet();
1778 :
1779 : // Set record format for calibration table application information
1780 0 : RecordDesc applyparDesc;
1781 0 : applyparDesc.addField ("t", TpDouble);
1782 0 : applyparDesc.addField ("table", TpString);
1783 0 : applyparDesc.addField ("interp", TpString);
1784 0 : applyparDesc.addField ("spw", TpArrayInt);
1785 0 : applyparDesc.addField ("field", TpArrayInt);
1786 0 : applyparDesc.addField ("calwt",TpBool);
1787 0 : applyparDesc.addField ("spwmap",TpArrayInt);
1788 0 : applyparDesc.addField ("opacity",TpFloat);
1789 :
1790 : // Create record with the requisite field values
1791 0 : Record applypar(applyparDesc);
1792 0 : applypar.define ("t", t);
1793 0 : applypar.define ("table", table);
1794 0 : applypar.define ("interp", interp);
1795 0 : applypar.define ("spw",Vector<Int>());
1796 0 : applypar.define ("field",Vector<Int>());
1797 : // applypar.define ("spw",getSpwIdx(spw));
1798 : // applypar.define ("field",getFieldIdx(field));
1799 0 : applypar.define ("calwt",calwt);
1800 0 : applypar.define ("spwmap",spwmap);
1801 0 : applypar.define ("opacity", opacity);
1802 :
1803 0 : String upType=type;
1804 0 : if (upType=="")
1805 : // Get type from table
1806 0 : upType = calTableType(table);
1807 :
1808 : // Must be upper case
1809 0 : upType.upcase();
1810 :
1811 : os << LogIO::NORMAL
1812 : << "Arranging to CORRUPT with:"
1813 0 : << LogIO::POST;
1814 :
1815 : // Add a new VisCal to the apply list
1816 0 : vc = createVisCal(upType,*vs_p);
1817 :
1818 0 : vc->setApply(applypar);
1819 :
1820 : os << LogIO::NORMAL << ". "
1821 0 : << vc->applyinfo()
1822 0 : << LogIO::POST;
1823 :
1824 0 : } catch (AipsError x) {
1825 : os << LogIO::SEVERE << x.getMesg()
1826 : << " Check inputs and try again."
1827 0 : << LogIO::POST;
1828 0 : if (vc) delete vc;
1829 0 : throw(AipsError("Error in Simulator::setapply."));
1830 : return false;
1831 : }
1832 :
1833 : // Creation apparently successful, so add to the apply list
1834 : // TBD: consolidate with above?
1835 : try {
1836 :
1837 0 : uInt napp=vc_p.nelements();
1838 0 : vc_p.resize(napp+1,false,true);
1839 0 : vc_p[napp] = vc;
1840 0 : vc=NULL;
1841 :
1842 : // Maintain sort of apply list
1843 0 : ve_p.setapply(vc_p);
1844 :
1845 0 : return true;
1846 :
1847 0 : } catch (AipsError x) {
1848 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1849 0 : << LogIO::POST;
1850 0 : if (vc) delete vc;
1851 0 : throw(AipsError("Error in Simulator::setapply."));
1852 : return false;
1853 : }
1854 : return false;
1855 : }
1856 :
1857 :
1858 :
1859 :
1860 :
1861 :
1862 : // Bool Simulator::corrupt() {
1863 16 : Bool Simulator::corrupt() {
1864 :
1865 : // VIS-plane (only) corruption
1866 :
1867 32 : LogIO os(LogOrigin("Simulator", "corrupt()", WHERE));
1868 :
1869 : try {
1870 :
1871 16 : ms_p->lock();
1872 16 : if(mssel_p) mssel_p->lock();
1873 :
1874 : // makeVisSet();
1875 : // AlwaysAssert(vs_p, AipsError);
1876 :
1877 : // Arrange the sort for efficient cal apply (different from sort order
1878 : // created by makeVisSet) and if there's a vs_p delete it and replace with
1879 : // this one. also delete it at the end of this routine
1880 32 : Block<Int> columns;
1881 : // include scan iteration
1882 16 : columns.resize(5);
1883 16 : columns[0]=MS::ARRAY_ID;
1884 16 : columns[1]=MS::SCAN_NUMBER;
1885 16 : columns[2]=MS::FIELD_ID;
1886 16 : columns[3]=MS::DATA_DESC_ID;
1887 16 : columns[4]=MS::TIME;
1888 :
1889 16 : if(vs_p) {
1890 16 : delete vs_p; vs_p=0;
1891 : }
1892 32 : Matrix<Int> noselection;
1893 16 : if(mssel_p) {
1894 16 : vs_p = new VisSet(*mssel_p,columns,noselection);
1895 : }
1896 : else {
1897 0 : vs_p = new VisSet(*ms_p,columns,noselection);
1898 : }
1899 16 : AlwaysAssert(vs_p, AipsError);
1900 :
1901 : // initCalSet() happens in VS creation unless there is a stored channel selection
1902 : // in the ms, and it equals the channel selection passed from here in mssel_p
1903 : // vs_p->initCalSet();
1904 :
1905 16 : VisIter& vi=vs_p->iter();
1906 32 : VisBuffer vb(vi);
1907 :
1908 : // Ensure VisEquation is ready - this sorts the VCs
1909 : // if we want a different order for corruption we will either need to
1910 : // implement the sort here or create a VE::setcorrupt(vc_p)
1911 16 : ve_p.setapply(vc_p);
1912 :
1913 : // set to corrupt Model down to B (T,D,G,etc) and correct Observed with AN,M,K
1914 16 : ve_p.setPivot(VisCal::B);
1915 :
1916 : // Apply
1917 16 : if (vc_p.nelements()>0) {
1918 : os << LogIO::NORMAL << "Doing visibility corruption."
1919 16 : << LogIO::POST;
1920 48 : for (uInt i=0;i<vc_p.nelements();i++) {
1921 : // os << vc_p[i]->longTypeName() << endl << vc_p[i]->siminfo() << endl <<
1922 : // "spwok = " << vc_p[i]->spwOK() << LogIO::POST;
1923 32 : os << vc_p[i]->siminfo() << "spwok = " << vc_p[i]->spwOK();
1924 32 : if (vc_p[i]->type() >= ve_p.pivot())
1925 16 : os << " in corrupt mode." << endl << LogIO::POST;
1926 : else
1927 16 : os << " in correct mode." << endl << LogIO::POST;
1928 : }
1929 34 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1930 : // Only if
1931 18 : if (ve_p.spwOK(vi.spectralWindow())) {
1932 4550 : for (vi.origin(); vi.more(); vi++) {
1933 :
1934 : // corrupt model to pivot, correct data up to pivot
1935 4532 : ve_p.collapseForSim(vb);
1936 :
1937 : // Deposit corrupted visibilities into DATA
1938 : // vi.setVis(vb.modelVisCube(), VisibilityIterator::Observed);
1939 4532 : vi.setVis(vb.visCube(), VisibilityIterator::Observed);
1940 : // for now, Also deposit in corrected
1941 : // (until newmmssimulator doesn't make corrected anymore)
1942 : // actually we should have this check if corrected is there,
1943 : // and if it is for some reason, copy data into it.
1944 : // RI TODO Sim::corrupt check for existence of Corrected
1945 4532 : vi.setVis(vb.visCube(), VisibilityIterator::Corrected);
1946 :
1947 : // RI TODO is this 100% right?
1948 4532 : vi.setWeightMat(vb.weightMat());
1949 :
1950 : }
1951 : }
1952 : else
1953 0 : cout << "Encountered data spw " << vi.spectralWindow() << " for which there is no (simulated) calibration." << endl;
1954 : }
1955 : }
1956 :
1957 :
1958 : // Old-fashioned noise, for now
1959 16 : if(ac_p != NULL){
1960 : // os << LogIO::WARN << "Using deprecated ACoh Noise - this will dissapear in the future - please switch to sm.setnoise" << LogIO::POST;
1961 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1962 0 : for (vi.origin(); vi.more(); vi++) {
1963 :
1964 : // affects vb.visibility() i.e. Observed
1965 0 : ac_p->apply(vb);
1966 0 : vi.setVis(vb.visibility(), VisibilityIterator::Observed);
1967 0 : vi.setVis(vb.visibility(), VisibilityIterator::Corrected);
1968 : }
1969 : }
1970 : }
1971 :
1972 : // Flush to disk
1973 16 : vs_p->flush();
1974 :
1975 : // since we made a specially sorted vs_p for corrupt, should we delete it
1976 : // and make one with the regular order?
1977 : // if(vs_p) {
1978 : // delete vs_p; vs_p=0;
1979 : // makeVisSet()
1980 : // }
1981 :
1982 16 : ms_p->relinquishAutoLocks();
1983 16 : ms_p->unlock();
1984 16 : if(mssel_p) mssel_p->unlock();
1985 :
1986 0 : } catch (std::exception& x) {
1987 0 : ms_p->relinquishAutoLocks();
1988 0 : ms_p->unlock();
1989 0 : if(mssel_p) mssel_p->unlock();
1990 0 : throw(AipsError(x.what()));
1991 : return false;
1992 : }
1993 16 : return true;
1994 : }
1995 :
1996 :
1997 :
1998 :
1999 :
2000 :
2001 :
2002 :
2003 :
2004 :
2005 :
2006 :
2007 :
2008 :
2009 :
2010 :
2011 :
2012 :
2013 :
2014 :
2015 2 : Bool Simulator::observe(const String& sourcename,
2016 : const String& spwname,
2017 : const Quantity& startTime,
2018 : const Quantity& stopTime,
2019 : const Bool add_observation,
2020 : const Bool state_sig,
2021 : const Bool state_ref,
2022 : const double& state_cal,
2023 : const double& state_load,
2024 : const unsigned int state_sub_scan,
2025 : const String& state_obs_mode,
2026 : const String& observername,
2027 : const String& projectname)
2028 : {
2029 6 : LogIO os(LogOrigin("Simulator", "observe()", WHERE));
2030 :
2031 :
2032 : try {
2033 :
2034 2 : if(!feedsHaveBeenSet && !feedsInitialized) {
2035 0 : os << "Feeds have not been set - using default " << feedMode_p << LogIO::WARN;
2036 0 : sim_p->initFeeds(feedMode_p);
2037 0 : feedsInitialized = true;
2038 : }
2039 2 : if(!timesHaveBeenSet_p) {
2040 : os << "Times have not been set - using defaults " << endl
2041 : << " Times will be interpreted as hour angles for first source"
2042 0 : << LogIO::WARN;
2043 : }
2044 :
2045 2 : sim_p->observe(sourcename, spwname, startTime, stopTime,
2046 : add_observation, state_sig, state_ref, state_cal,state_load,state_sub_scan,state_obs_mode,observername,projectname);
2047 :
2048 :
2049 2 : if(ms_p) delete ms_p; ms_p=0;
2050 2 : if(mssel_p) delete mssel_p; mssel_p=0;
2051 4 : ms_p = new MeasurementSet(msname_p,
2052 2 : TableLock(TableLock::AutoNoReadLocking),
2053 2 : Table::Update);
2054 :
2055 2 : ms_p->flush();
2056 2 : ms_p->unlock();
2057 :
2058 0 : } catch (AipsError x) {
2059 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2060 0 : return false;
2061 : }
2062 2 : return true;
2063 : }
2064 :
2065 :
2066 :
2067 :
2068 27 : Bool Simulator::observemany(const Vector<String>& sourcenames,
2069 : const String& spwname,
2070 : const Vector<Quantity>& startTimes,
2071 : const Vector<Quantity>& stopTimes,
2072 : const Vector<MDirection>& directions,
2073 : const Bool add_observation=true,
2074 : const Bool state_sig=true,
2075 : const Bool state_ref=true,
2076 : const double& state_cal=0.,
2077 : const double& state_load=0.,
2078 : const unsigned int state_sub_scan=0,
2079 : const String& state_obs_mode="OBSERVE_TARGET#ON_SOURCE",
2080 : const String& observername="CASA simulator",
2081 : const String& projectname="CASA simulation")
2082 : {
2083 81 : LogIO os(LogOrigin("Simulator", "observemany()", WHERE));
2084 :
2085 :
2086 : try {
2087 :
2088 27 : if(!feedsHaveBeenSet && !feedsInitialized) {
2089 0 : os << "Feeds have not been set - using default " << feedMode_p << LogIO::WARN;
2090 0 : sim_p->initFeeds(feedMode_p);
2091 0 : feedsInitialized = true;
2092 : }
2093 27 : if(!timesHaveBeenSet_p) {
2094 : os << "Times have not been set - using defaults " << endl
2095 : << " Times will be interpreted as hour angles for first source"
2096 0 : << LogIO::WARN;
2097 : }
2098 :
2099 27 : sim_p->observe(sourcenames, spwname, startTimes, stopTimes, directions,
2100 : add_observation, state_sig, state_ref, state_cal,state_load,state_sub_scan,state_obs_mode,observername,projectname);
2101 :
2102 27 : if(ms_p) delete ms_p; ms_p=0;
2103 27 : if(mssel_p) delete mssel_p; mssel_p=0;
2104 54 : ms_p = new MeasurementSet(msname_p,
2105 27 : TableLock(TableLock::AutoNoReadLocking),
2106 27 : Table::Update);
2107 :
2108 27 : ms_p->flush();
2109 27 : ms_p->unlock();
2110 :
2111 0 : } catch (AipsError x) {
2112 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2113 0 : return false;
2114 : }
2115 27 : return true;
2116 : }
2117 :
2118 :
2119 :
2120 :
2121 :
2122 :
2123 27 : Bool Simulator::predict(const Vector<String>& modelImage,
2124 : const String& compList,
2125 : const Bool incremental) {
2126 :
2127 81 : LogIO os(LogOrigin("Simulator", "predict()", WHERE));
2128 :
2129 : // Note that incremental here does not apply to se_p->predict(false),
2130 : // Rather it means: add the calculated model visibility to the data visibility.
2131 : // We return a MS with Data, Model, and Corrected columns identical
2132 :
2133 : try {
2134 :
2135 : os << "Predicting visibilities using model: " << modelImage <<
2136 27 : " and componentList: " << compList << LogIO::POST;
2137 27 : if (incremental) {
2138 0 : os << "The data column will be incremented" << LogIO::POST;
2139 : } else {
2140 27 : os << "The data column will be replaced" << LogIO::POST;
2141 : }
2142 27 : if(!ms_p) {
2143 : os << "MeasurementSet pointer is null : logic problem!"
2144 0 : << LogIO::EXCEPTION;
2145 : }
2146 :
2147 27 : bool heterogeneous=False;
2148 744 : for (uInt i=1;i<diam_p.nelements();i++)
2149 717 : if (diam_p(i)!=diam_p(0)) heterogeneous=True;
2150 27 : if (heterogeneous) {
2151 0 : if (compList!=String("")) {
2152 0 : os<<"Heterogeneous array is not supported for simulation from components. Primary beam will be applied by telescope name.\n"<<LogIO::WARN;
2153 : }
2154 0 : if ((modelImage[0]!=String(""))&&(ftmachine_p!="mosaic")) {
2155 0 : os<<"Heterogeneous array is only supported for ALMA,ACA,OVRO telescopes, and only with option ftmachine=mosaic."<<LogIO::WARN;
2156 : }
2157 : }
2158 :
2159 27 : ms_p->lock();
2160 27 : if(mssel_p) mssel_p->lock();
2161 27 : if (!createSkyEquation( modelImage, compList)) {
2162 0 : os << LogIO::SEVERE << "Failed to create SkyEquation" << LogIO::POST;
2163 0 : return false;
2164 : }
2165 27 : if (incremental) {
2166 0 : se_p->predict(false,MS::MODEL_DATA);
2167 : } else {
2168 27 : se_p->predict(false,MS::DATA); //20091030 RI changed SE::predict to use DATA
2169 : }
2170 27 : destroySkyEquation();
2171 :
2172 : // Copy the predicted visibilities over to the observed and
2173 : // the corrected data columns
2174 27 : makeVisSet();
2175 :
2176 27 : VisIter& vi = vs_p->iter();
2177 54 : VisBuffer vb(vi);
2178 27 : vi.origin();
2179 27 : vi.originChunks();
2180 :
2181 : //os << "Copying predicted visibilities from MODEL_DATA to DATA and CORRECTED_DATA" << LogIO::POST;
2182 :
2183 406 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()){
2184 1008 : for (vi.origin(); vi.more(); vi++) {
2185 : // vb.setVisCube(vb.modelVisCube());
2186 :
2187 629 : if (incremental) {
2188 0 : vi.setVis( (vb.modelVisCube() + vb.visCube()),
2189 0 : VisibilityIterator::Corrected);
2190 0 : vi.setVis(vb.correctedVisCube(),VisibilityIterator::Observed);
2191 :
2192 : // model=1 is more consistent with VS::initCalSet
2193 : // vi.setVis(vb.correctedVisCube(),VisibilityIterator::Model);
2194 : } else {
2195 : // from above, the prediction is now already in Observed.
2196 : // RI TODO remove scratch columns from NewMSSimulator;
2197 : // until then we;ll just leave them 1 and Corr=Obs (for imaging)
2198 : //vi.setVis(vb.visCube(),VisibilityIterator::Observed);
2199 629 : vi.setVis(vb.visCube(),VisibilityIterator::Corrected);
2200 : }
2201 629 : vb.setModelVisCube(Complex(1.0,0.0));
2202 629 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
2203 : }
2204 : }
2205 27 : ms_p->unlock();
2206 27 : if(mssel_p) mssel_p->lock();
2207 :
2208 0 : } catch (AipsError x) {
2209 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2210 0 : ms_p->unlock();
2211 0 : if(mssel_p) mssel_p->lock();
2212 0 : return false;
2213 : }
2214 27 : return true;
2215 : }
2216 :
2217 :
2218 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2219 : // copied from SynthesisImager
2220 27 : void Simulator::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2221 : {
2222 81 : LogIO os(LogOrigin("Simulator", "getVPRecord",WHERE));
2223 :
2224 27 : VPManager *vpman=VPManager::Instance();
2225 27 : if ((doDefaultVP_p)||(vpTableStr_p != String(""))) {
2226 10 : if (doDefaultVP_p) {
2227 10 : os << "Using default Voltage Patterns from the VPManager" << LogIO::POST;
2228 10 : vpman->reset();
2229 : } else {
2230 0 : os << "Loading Voltage Pattern information from " << vpTableStr_p << LogIO::POST;
2231 0 : vpman->loadfromtable(vpTableStr_p );
2232 : }
2233 10 : os << "Temporary alert : The state of the vpmanager tool has been modified by loading these primary beam models. If any of your scripts rely on the vpmanager state being preserved throughout your CASA session, please use vp.saveastable() and vp.loadfromtable() as needed. "<< LogIO::POST;
2234 : }
2235 : // if dodefault=false, and no table is provided, it will used what is in
2236 : // the vpmanager already.
2237 : else {
2238 17 : os << "Using Voltage Patterns from the VPManager" << LogIO::POST;
2239 : }
2240 :
2241 : // PBMath::CommonPB kpb;
2242 27 : PBMath::enumerateCommonPB(telescop, kpb);
2243 : // Record rec;
2244 27 : vpman->getvp(rec, telescop);
2245 :
2246 : //ostringstream oss; rec.print(oss);
2247 : //os << "Using VP model : " << oss.str() << LogIO::POST;
2248 :
2249 27 : if(rec.empty()){
2250 0 : if((telescop=="EVLA")){
2251 0 : os << LogIO::WARN << "vpmanager does not list EVLA. Using VLA beam parameters. To use EVLA beams, please use the vpmanager and set the vptable parameter in tclean (see inline help for vptable)." << LogIO::POST;
2252 0 : telescop="VLA";
2253 0 : vpman->getvp(rec, telescop);
2254 0 : kpb=PBMath::VLA;
2255 : }
2256 : else{
2257 0 : os << LogIO::WARN << "vpmanager does not have a beam for "+telescop <<".\n Please use the vpmanager to define one (and optionally save its state as a table and supply its name via the vptable parameter.)" << LogIO::POST;
2258 0 : kpb=PBMath::UNKNOWN;
2259 : }
2260 : }
2261 :
2262 27 : }// get VPRecord
2263 :
2264 :
2265 :
2266 27 : Bool Simulator::createSkyEquation(const Vector<String>& image,
2267 : const String complist)
2268 : {
2269 :
2270 81 : LogIO os(LogOrigin("Simulator", "createSkyEquation()", WHERE));
2271 :
2272 : try {
2273 27 : if(sm_p==0) {
2274 27 : sm_p = new CleanImageSkyModel();
2275 : }
2276 27 : AlwaysAssert(sm_p, AipsError);
2277 :
2278 : // Add the componentlist
2279 27 : if(complist!="") {
2280 7 : if(!Table::isReadable(complist)) {
2281 : os << LogIO::SEVERE << "ComponentList " << complist
2282 0 : << " not readable" << LogIO::POST;
2283 0 : return false;
2284 : }
2285 7 : componentList_p=new ComponentList(complist, true);
2286 7 : if(componentList_p==0) {
2287 : os << LogIO::SEVERE << "Cannot create ComponentList from " << complist
2288 0 : << LogIO::POST;
2289 0 : return false;
2290 : }
2291 7 : if(!sm_p->add(*componentList_p)) {
2292 : os << LogIO::SEVERE << "Cannot add ComponentList " << complist
2293 0 : << " to SkyModel" << LogIO::POST;
2294 0 : return false;
2295 : }
2296 : } else {
2297 20 : componentList_p=0;
2298 : }
2299 :
2300 27 : nmodels_p = image.nelements();
2301 27 : if (nmodels_p == 1 && image(0) == "") nmodels_p = 0;
2302 :
2303 27 : int models_found=0;
2304 27 : if (nmodels_p > 0) {
2305 23 : images_p.resize(nmodels_p);
2306 :
2307 46 : for (Int model=0;model<Int(nmodels_p);model++) {
2308 23 : if(image(model)=="") {
2309 : os << LogIO::SEVERE << "Need a name for model " << model+1
2310 0 : << LogIO::POST;
2311 0 : return false;
2312 : } else {
2313 23 : if(!Table::isReadable(image(model))) {
2314 0 : os << LogIO::SEVERE << image(model) << " is unreadable" << LogIO::POST;
2315 : } else {
2316 23 : images_p[model]=0;
2317 : os << "Opening model " << model << " named "
2318 23 : << image(model) << LogIO::POST;
2319 23 : images_p[model]=new PagedImage<Float>(image(model));
2320 23 : AlwaysAssert(images_p[model], AipsError);
2321 : // RI TODO is this a logic problem with more than one source??
2322 : // Add distance
2323 23 : if((distance_p.nelements() > 0 && distance_p.nelements() <= static_cast<uInt>(nField)) && abs(distance_p[nField-1].get().getValue())>0.0) {
2324 0 : os << " Refocusing to distance " << distance_p[nField-1].get("km").getValue()
2325 0 : << " km" << LogIO::POST;
2326 0 : Record info(images_p[model]->miscInfo());
2327 0 : info.define("distance", distance_p[nField-1].get("m").getValue());
2328 0 : images_p[model]->setMiscInfo(info);
2329 : }
2330 :
2331 : // FTMachine only works in Hz and LSRK
2332 23 : CoordinateSystem cs = images_p[model]->coordinates();
2333 23 : String errorMsg;
2334 23 : cs.setSpectralConversion(errorMsg,MFrequency::showType(MFrequency::LSRK));
2335 23 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
2336 23 : AlwaysAssert(spectralIndex>=0, AipsError);
2337 23 : SpectralCoordinate spectralCoord=cs.spectralCoordinate(spectralIndex);
2338 23 : Vector<String> units(1); units = "Hz";
2339 : // doesn't work: cs.spectralCoordinate(spectralIndex).setWorldAxisUnits(units);
2340 23 : spectralCoord.setWorldAxisUnits(units);
2341 : // put spectralCoord back into cs
2342 23 : cs.replaceCoordinate(spectralCoord,spectralIndex);
2343 : // and cs into image
2344 23 : images_p[model]->setCoordinateInfo(cs);
2345 :
2346 :
2347 23 : if(sm_p->add(*images_p[model])!=model) {
2348 0 : os << LogIO::SEVERE << "Error adding model " << model << LogIO::POST;
2349 0 : return false;
2350 : }
2351 23 : models_found++;
2352 : }
2353 : }
2354 : }
2355 : }
2356 :
2357 27 : if(models_found<=0 and componentList_p==0) {
2358 0 : os << LogIO::SEVERE << "No model images found" << LogIO::POST;
2359 0 : return false;
2360 : }
2361 :
2362 :
2363 :
2364 27 : if(vs_p) {
2365 27 : delete vs_p; vs_p=0;
2366 : }
2367 27 : makeVisSet();
2368 :
2369 27 : cft_p = new SimpleComponentFTMachine();
2370 :
2371 27 : MeasurementSet *ams=0;
2372 :
2373 27 : if(mssel_p) {
2374 27 : ams=mssel_p;
2375 : }
2376 : else {
2377 0 : ams=ms_p;
2378 : }
2379 : // 2016 from SynthesisImager:
2380 54 : MSColumns msc(*ams);
2381 54 : String telescop=msc.observation().telescopeName()(0);
2382 : PBMath::CommonPB kpb;
2383 54 : Record rec;
2384 27 : getVPRecord( rec, kpb, telescop );
2385 :
2386 27 : if((ftmachine_p=="sd")||(ftmachine_p=="both")||(ftmachine_p=="mosaic")) {
2387 10 : if(!gvp_p) {
2388 10 : if (rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2389 10 : String commonPBName;
2390 10 : PBMath::nameCommonPB(kpb, commonPBName);
2391 10 : os << "Using common PB "<<commonPBName<<" for beam calculation for telescope " << telescop << LogIO::POST;
2392 10 : gvp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p, skyPosThreshold_p);
2393 : }
2394 : else{
2395 0 : PBMath myPB(rec);
2396 0 : String whichPBMath;
2397 0 : PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2398 0 : os << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2399 0 : gvp_p= new VPSkyJones(telescop, myPB, parAngleInc_p, squintType_p, skyPosThreshold_p);
2400 0 : kpb=PBMath::DEFAULT;
2401 : }
2402 : }
2403 10 : if(ftmachine_p=="sd") {
2404 10 : os << "Performing Single dish gridding " << LogIO::POST;
2405 10 : if(gridfunction_p=="pb") {
2406 10 : ft_p = new SDGrid(*gvp_p, cache_p/2, tile_p, gridfunction_p);
2407 : }
2408 : else {
2409 0 : ft_p = new SDGrid(cache_p/2, tile_p, gridfunction_p);
2410 : }
2411 : }
2412 0 : else if(ftmachine_p=="mosaic") {
2413 : // RI TODO need stokesString for current spw - e.g. currSpw()?
2414 0 : os << "Performing Mosaic gridding for model images (uv convolution)" << LogIO::POST;
2415 : //2016 from SynthesisImager:
2416 0 : ft_p = new MosaicFTNew(gvp_p, mLocation_p, stokesString_p[0], cache_p/2, tile_p, true);
2417 0 : PBMathInterface::PBClass pbtype=PBMathInterface::AIRY;
2418 0 : if(rec.asString("name")=="IMAGE")
2419 0 : pbtype=PBMathInterface::IMAGE;
2420 : ///Use Heterogenous array mode for the following
2421 0 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
2422 0 : || (kpb==PBMath::ALMA)){
2423 0 : String PBName;
2424 0 : PBMath::nameCommonPB(kpb,PBName);
2425 0 : os << "Enabling Heterogeneous Array for PBMath "<<PBName<<LogIO::POST;
2426 : // Will use Airys - to do something more complex we need to
2427 : // use HetArrayConvFunv::fromRecord
2428 0 : if ((kpb==PBMath::ACA) || (kpb==PBMath::ALMA)) {
2429 0 : os << "For model images, will use 10.7m Airy PB for diameter=12 dishes, and 6.25m Airy PB for diameter=7 dishes." << LogIO::POST;
2430 : } else {
2431 0 : os << "For model images, will use Airy PB scaled to dish diameter."<<LogIO::POST;
2432 : }
2433 0 : CountedPtr<SimplePBConvFunc> mospb=new HetArrayConvFunc(pbtype, "");
2434 0 : static_cast<MosaicFTNew &>(*ft_p).setConvFunc(mospb);
2435 : }
2436 : //gvp_p->summary();
2437 : ///2016
2438 : }
2439 0 : else if(ftmachine_p=="both") {
2440 : os << "Performing single dish gridding with convolution function "
2441 0 : << gridfunction_p << LogIO::POST;
2442 : os << "and interferometric gridding with convolution function SF"
2443 0 : << LogIO::POST;
2444 :
2445 0 : ft_p = new GridBoth(*gvp_p, cache_p/2, tile_p,
2446 0 : mLocation_p,
2447 0 : gridfunction_p, "SF", padding_p);
2448 : }
2449 :
2450 10 : VisIter& vi(vs_p->iter());
2451 : // Get bigger chunks o'data: this should be tuned some time
2452 : // since it may be wrong for e.g. spectral line
2453 10 : vi.setRowBlocking(100);
2454 : }
2455 :
2456 : else {
2457 17 : os << "Synthesis gridding " << LogIO::POST;
2458 : // Now make the FTMachine
2459 : // if(wprojPlanes_p>1) {
2460 17 : if (ftmachine_p=="wproject") {
2461 0 : os << "Fourier transforms will use specified common tangent point:" << LogIO::POST;
2462 : // RI TODO how does this work with more than one field?
2463 0 : os << formatDirection(sourceDirection_p[nField-1]) << LogIO::POST;
2464 : // ft_p = new WProjectFT(*ams, facets_p, cache_p/2, tile_p, false);
2465 0 : ft_p = new WProjectFT(wprojPlanes_p, mLocation_p, cache_p/2, tile_p, false);
2466 : }
2467 17 : else if (ftmachine_p=="pbwproject") {
2468 : os << "Fourier transforms will use specified common tangent point and PBs"
2469 0 : << LogIO::POST;
2470 0 : os << formatDirection(sourceDirection_p[nField-1]) << LogIO::POST;
2471 :
2472 : // if (!epJ_p)
2473 : os << "Antenna pointing related term (EPJones) not set. "
2474 : << "This is required when using pbwproject FTMachine."
2475 : << "(gmoellen 06Nov20: pointing errors temporarily disabled)"
2476 0 : << LogIO::EXCEPTION;
2477 :
2478 : /*
2479 : doVP_p = false; // Since this FTMachine includes PB
2480 : if (wprojPlanes_p<=1)
2481 : {
2482 : os << LogIO::NORMAL
2483 : << "You are using wprojplanes=1. Doing co-planar imaging "
2484 : << "(no w-projection needed)"
2485 : << LogIO::POST;
2486 : os << "Performing pb-projection"
2487 : << LogIO::POST;
2488 : }
2489 : if((wprojPlanes_p>1)&&(wprojPlanes_p<64))
2490 : {
2491 : os << LogIO::WARN
2492 : << "No. of w-planes set too low for W projection - recommend at least 128"
2493 : << LogIO::POST;
2494 : os << "Performing pb + w-plane projection"
2495 : << LogIO::POST;
2496 : }
2497 : // epJ_p = new EPJones(*vs_p);
2498 : // epJ_p->load(epJTableName_p,"","diagonal");
2499 : if(!gvp_p)
2500 : {
2501 : os << "Using defaults for primary beams used in gridding" << LogIO::POST;
2502 : gvp_p=new VPSkyJones(*ms_p, true, parAngleInc_p, squintType_p);
2503 : }
2504 : // ft_p = new PBWProjectFT(*ms_p, epJ, gvp_p, facets_p, cache_p/2,
2505 : // doPointing, tile_p, paStep_p,
2506 : // pbLimit_p, true);
2507 :
2508 : String cfCacheDirName = "cache";
2509 : if (mssel_p)
2510 : ft_p = new PBWProjectFT(*mssel_p, epJ_p,
2511 : // gvp_p,
2512 : wprojPlanes_p, cache_p/2,
2513 : cfCacheDirName,
2514 : applyPointingOffsets_p, doPBCorrection_p,
2515 : tile_p,
2516 : 0.0, // Not required here. parAngleInc_p is used in gvp_p
2517 : pbLimit_p, true);
2518 : else
2519 : ft_p = new PBWProjectFT(*ms_p, epJ_p,
2520 : // gvp_p,
2521 : wprojPlanes_p, cache_p/2,
2522 : cfCacheDirName,
2523 : applyPointingOffsets_p, doPBCorrection_p,
2524 : tile_p,
2525 : 0.0, // Not required here. parAngleInc_p is used in gvp_p
2526 : pbLimit_p, true);
2527 : AlwaysAssert(ft_p, AipsError);
2528 : cft_p = new SimpleComponentFTMachine();
2529 : AlwaysAssert(cft_p, AipsError);
2530 :
2531 : */
2532 : }
2533 : else {
2534 17 : os << "Fourier transforms will use image centers as tangent points" << LogIO::POST;
2535 17 : ft_p = new GridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p, padding_p);
2536 : }
2537 : }
2538 27 : AlwaysAssert(ft_p, AipsError);
2539 :
2540 : // tell ftmachine about the transformations in model images above - no.
2541 : //Vector<Int> dataspectralwindowids_p;
2542 : //Bool freqFrameValid_p = true;
2543 : //ft_p->setSpw(dataspectralwindowids_p, freqFrameValid_p);
2544 :
2545 27 : se_p = new SkyEquation ( *sm_p, *vs_p, *ft_p, *cft_p );
2546 :
2547 : // Now add any SkyJones that are needed
2548 : // (vp_p is not applied to images if ftmachine=mosaic)
2549 27 : if(doVP_p) {
2550 17 : if (rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2551 : // String commonPBName;
2552 : // PBMath::nameCommonPB(kpb, commonPBName);
2553 : // os << "SkyEquation using common PB "<<commonPBName<<" for beam calculation for telescope " << telescop << LogIO::POST;
2554 17 : vp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p, skyPosThreshold_p);
2555 : } else {
2556 0 : PBMath myPB(rec);
2557 : // String whichPBMath;
2558 : // PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2559 : // os << "SkyEquation using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2560 0 : vp_p= new VPSkyJones(telescop, myPB, parAngleInc_p, squintType_p, skyPosThreshold_p);
2561 0 : kpb=PBMath::DEFAULT;
2562 : }
2563 17 : vp_p->summary();
2564 17 : se_p->setSkyJones(*vp_p);
2565 : } else {
2566 10 : vp_p=0;
2567 : }
2568 0 : } catch (AipsError x) {
2569 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2570 0 : ms_p->unlock();
2571 0 : if(mssel_p) mssel_p->lock();
2572 0 : return false;
2573 : }
2574 27 : return true;
2575 : };
2576 :
2577 27 : void Simulator::destroySkyEquation()
2578 : {
2579 27 : if(se_p) delete se_p; se_p=0;
2580 27 : if(sm_p) delete sm_p; sm_p=0;
2581 27 : if(vp_p) delete vp_p; vp_p=0;
2582 27 : if(componentList_p) delete componentList_p; componentList_p=0;
2583 :
2584 50 : for (Int model=0;model<Int(nmodels_p);model++) {
2585 23 : if(images_p[model]) delete images_p[model]; images_p[model]=0;
2586 : }
2587 27 : };
2588 :
2589 :
2590 :
2591 :
2592 :
2593 :
2594 :
2595 :
2596 :
2597 :
2598 :
2599 :
2600 :
2601 :
2602 0 : Bool Simulator::setlimits(const Double shadowLimit,
2603 : const Quantity& elevationLimit)
2604 : {
2605 :
2606 0 : LogIO os(LogOrigin("Simulator", "setlimits()", WHERE));
2607 :
2608 : try {
2609 :
2610 0 : sim_p->setFractionBlockageLimit( shadowLimit );
2611 0 : sim_p->setElevationLimit( elevationLimit );
2612 0 : } catch (AipsError x) {
2613 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2614 0 : return false;
2615 : }
2616 0 : return true;
2617 : }
2618 :
2619 28 : Bool Simulator::setauto(const Double autocorrwt)
2620 : {
2621 :
2622 56 : LogIO os(LogOrigin("Simulator", "setauto()", WHERE));
2623 :
2624 : try {
2625 :
2626 28 : sim_p->setAutoCorrelationWt(autocorrwt);
2627 :
2628 : } catch (AipsError x) {
2629 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2630 : return false;
2631 : }
2632 28 : return true;
2633 : }
2634 :
2635 129 : void Simulator::makeVisSet() {
2636 :
2637 129 : if(vs_p) return;
2638 :
2639 140 : Block<int> sort(0);
2640 70 : sort.resize(5);
2641 70 : sort[0] = MS::FIELD_ID;
2642 70 : sort[1] = MS::FEED1;
2643 70 : sort[2] = MS::ARRAY_ID;
2644 70 : sort[3] = MS::DATA_DESC_ID;
2645 70 : sort[4] = MS::TIME;
2646 140 : Matrix<Int> noselection;
2647 70 : if(mssel_p) {
2648 70 : vs_p = new VisSet(*mssel_p,sort,noselection);
2649 : }
2650 : else {
2651 0 : if (ms_p)
2652 0 : vs_p = new VisSet(*ms_p,sort,noselection);
2653 : else
2654 0 : throw(AipsError("No measurement set open in Simulator."));
2655 : }
2656 70 : AlwaysAssert(vs_p, AipsError);
2657 :
2658 : }
2659 :
2660 :
2661 :
2662 10 : Bool Simulator::setoptions(const String& ftmachine, const Int cache, const Int tile,
2663 : const String& gridfunction, const MPosition& mLocation,
2664 : const Float padding, const Int facets, const Double maxData,
2665 : const Int wprojPlanes)
2666 : {
2667 20 : LogIO os(LogOrigin("Simulator", "setoptions()", WHERE));
2668 :
2669 10 : os << "Setting processing options" << LogIO::POST;
2670 :
2671 10 : sim_p->setMaxData(maxData*1024.0*1024.0);
2672 :
2673 10 : ftmachine_p=downcase(ftmachine);
2674 10 : if(cache>0) cache_p=cache;
2675 10 : if(tile>0) tile_p=tile;
2676 10 : gridfunction_p=downcase(gridfunction);
2677 10 : mLocation_p=mLocation;
2678 10 : if(padding>=1.0) {
2679 10 : padding_p=padding;
2680 : }
2681 10 : facets_p=facets;
2682 10 : wprojPlanes_p = wprojPlanes;
2683 : // Destroy the FTMachine
2684 10 : if(ft_p) {delete ft_p; ft_p=0;}
2685 10 : if(cft_p) {delete cft_p; cft_p=0;}
2686 :
2687 20 : return true;
2688 : }
2689 :
2690 :
2691 0 : Bool Simulator::detached() const
2692 : {
2693 0 : return false;
2694 : }
2695 :
2696 0 : String Simulator::formatDirection(const MDirection& direction) {
2697 0 : MVAngle mvRa=direction.getAngle().getValue()(0);
2698 0 : MVAngle mvDec=direction.getAngle().getValue()(1);
2699 0 : ostringstream oss;
2700 0 : oss.setf(ios::left, ios::adjustfield);
2701 0 : oss.width(14);
2702 0 : oss << mvRa(0.0).string(MVAngle::TIME,8);
2703 0 : oss.width(14);
2704 0 : oss << mvDec.string(MVAngle::DIG2,8);
2705 0 : oss << " " << MDirection::showType(direction.getRefPtr()->getType());
2706 0 : return String(oss);
2707 : }
2708 :
2709 0 : String Simulator::formatTime(const Double time) {
2710 0 : MVTime mvtime(Quantity(time, "s"));
2711 0 : return mvtime.string(MVTime::DMY,7);
2712 : }
2713 :
2714 :
2715 :
2716 43 : Bool Simulator::setdata(const Vector<Int>& spectralwindowids,
2717 : const Vector<Int>& fieldids,
2718 : const String& msSelect)
2719 :
2720 : {
2721 :
2722 :
2723 129 : LogIO os(LogOrigin("Simulator", "setdata()", WHERE));
2724 :
2725 43 : if(!ms_p) {
2726 : os << LogIO::SEVERE << "Program logic error: MeasurementSet pointer ms_p not yet set"
2727 0 : << LogIO::POST;
2728 0 : return false;
2729 : }
2730 :
2731 : try {
2732 :
2733 43 : os << "Selecting data" << LogIO::POST;
2734 :
2735 : // Map the selected spectral window ids to data description ids
2736 86 : MSDataDescColumns dataDescCol(ms_p->dataDescription());
2737 86 : Vector<Int> ddSpwIds=dataDescCol.spectralWindowId().getColumn();
2738 :
2739 43 : Vector<Int> datadescids(0);
2740 86 : for (uInt row=0; row<ddSpwIds.nelements(); row++) {
2741 43 : Bool found=false;
2742 86 : for (uInt j=0; j<spectralwindowids.nelements(); j++) {
2743 43 : if (ddSpwIds(row)==spectralwindowids(j)) found=true;
2744 : };
2745 43 : if (found) {
2746 43 : datadescids.resize(datadescids.nelements()+1,true);
2747 43 : datadescids(datadescids.nelements()-1)=row;
2748 : };
2749 : };
2750 :
2751 43 : if(vs_p) delete vs_p; vs_p=0;
2752 43 : if(mssel_p) delete mssel_p; mssel_p=0;
2753 :
2754 : // If a selection has been made then close the current MS
2755 : // and attach to a new selected MS. We do this on the original
2756 : // MS.
2757 43 : if(fieldids.nelements()>0||datadescids.nelements()>0) {
2758 43 : os << "Performing selection on MeasurementSet" << LogIO::POST;
2759 43 : Table& original=*ms_p;
2760 :
2761 : // Now we make a condition to do the old FIELD_ID, SPECTRAL_WINDOW_ID
2762 : // selection
2763 86 : TableExprNode condition;
2764 86 : String colf=MS::columnName(MS::FIELD_ID);
2765 86 : String cols=MS::columnName(MS::DATA_DESC_ID);
2766 43 : if(fieldids.nelements()>0&&datadescids.nelements()>0){
2767 27 : condition=original.col(colf).in(fieldids)&&original.col(cols).in(datadescids);
2768 27 : os << "Selecting on field and spectral window ids" << LogIO::POST;
2769 : }
2770 16 : else if(datadescids.nelements()>0) {
2771 16 : condition=original.col(cols).in(datadescids);
2772 16 : os << "Selecting on spectral window id" << LogIO::POST;
2773 : }
2774 0 : else if(fieldids.nelements()>0) {
2775 0 : condition=original.col(colf).in(fieldids);
2776 0 : os << "Selecting on field id" << LogIO::POST;
2777 : }
2778 :
2779 : // Now remake the original ms
2780 43 : mssel_p = new MeasurementSet(original(condition));
2781 :
2782 : //AlwaysAssert(mssel_p, AipsError);
2783 : //mssel_p->rename(msname_p+"/SELECTED_TABLE", Table::Scratch);
2784 43 : if(mssel_p->nrow()==0) {
2785 0 : delete mssel_p; mssel_p=0;
2786 : os << LogIO::WARN
2787 : << "Selection is empty: reverting to original MeasurementSet"
2788 0 : << LogIO::POST;
2789 0 : mssel_p=new MeasurementSet(original);
2790 : }
2791 : else {
2792 43 : mssel_p->flush();
2793 : }
2794 :
2795 : }
2796 : else {
2797 0 : mssel_p=new MeasurementSet(*ms_p);
2798 : }
2799 : {
2800 43 : Int len = msSelect.length();
2801 43 : Int nspace = msSelect.freq (' ');
2802 43 : Bool nullSelect=(msSelect.empty() || nspace==len);
2803 43 : if (!nullSelect) {
2804 0 : os << "Now applying selection string " << msSelect << LogIO::POST;
2805 : MeasurementSet* mssel_p2;
2806 : // Apply the TAQL selection string, to remake the original MS
2807 0 : String parseString="select from $1 where " + msSelect;
2808 0 : mssel_p2=new MeasurementSet(tableCommand(parseString,*mssel_p).table());
2809 0 : AlwaysAssert(mssel_p2, AipsError);
2810 : // Rename the selected MS as */SELECTED_TABLE2
2811 : //mssel_p2->rename(msname_p+"/SELECTED_TABLE2", Table::Scratch);
2812 0 : if (mssel_p2->nrow()==0) {
2813 : os << LogIO::WARN
2814 : << "Selection string results in empty MS: "
2815 : << "reverting to original MeasurementSet"
2816 0 : << LogIO::POST;
2817 0 : delete mssel_p2;
2818 : } else {
2819 0 : if (mssel_p) {
2820 0 : delete mssel_p;
2821 0 : mssel_p=mssel_p2;
2822 0 : mssel_p->flush();
2823 : }
2824 : }
2825 : } else {
2826 43 : os << "No selection string given" << LogIO::POST;
2827 : }
2828 :
2829 43 : if(mssel_p->nrow()!=ms_p->nrow()) {
2830 1 : os << "By selection " << ms_p->nrow() << " rows are reduced to "
2831 1 : << mssel_p->nrow() << LogIO::POST;
2832 : }
2833 : else {
2834 42 : os << "Selection did not drop any rows" << LogIO::POST;
2835 : }
2836 : }
2837 :
2838 : // Now create the VisSet
2839 43 : if(vs_p) delete vs_p; vs_p=0;
2840 43 : makeVisSet();
2841 : //Now assign the source directions to something selected or sensible
2842 : {
2843 : //Int fieldsel=0;
2844 43 : if(fieldids.nelements() >0) {
2845 : //fieldsel=fieldids(0);
2846 : // RI TODO does sim:setdata need this?
2847 27 : nField=fieldids.nelements();
2848 514 : for (Int i=0;i<nField;i++) {
2849 : // RI TODO check whether index in field column is just i or need
2850 : // to search for fieldid
2851 487 : (vs_p->iter()).msColumns().field().name().get(i,sourceName_p[i]);
2852 487 : sourceDirection_p[i]=(vs_p->iter()).msColumns().field().phaseDirMeas(i);
2853 487 : (vs_p->iter()).msColumns().field().code().get(i,calCode_p[i]);
2854 : }
2855 : }
2856 : }
2857 43 : return true;
2858 0 : } catch (AipsError x) {
2859 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
2860 0 : << LogIO::POST;
2861 0 : return false;
2862 : }
2863 : return true;
2864 : }
2865 :
2866 : } // end namespace casa
|
