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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: 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 0 : msname_p(String("")), ms_p(0), mssel_p(0), vs_p(0),
106 0 : seed_p(11111),
107 0 : ac_p(0), distance_p(0), ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
108 0 : sim_p(0),
109 : // epJ_p(0),
110 0 : epJTableName_p(),
111 0 : prtlev_(0)
112 : {
113 0 : }
114 :
115 :
116 31 : Simulator::Simulator(String& msname)
117 31 : : msname_p(msname), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
118 31 : ac_p(0), distance_p(0),ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
119 31 : sim_p(0),
120 : // epJ_p(0),
121 31 : epJTableName_p(),
122 62 : prtlev_(0)
123 : {
124 31 : defaults();
125 :
126 31 : if(!sim_p) {
127 31 : sim_p= new NewMSSimulator(msname);
128 : //sim_p->setPrtlev(prtlev());
129 : }
130 :
131 31 : ve_p.setPrtlev(prtlev());
132 :
133 : // Make a MeasurementSet object for the disk-base MeasurementSet that we just
134 : // created
135 31 : ms_p = new MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
136 31 : Table::Update);
137 31 : AlwaysAssert(ms_p, AipsError);
138 :
139 31 : }
140 :
141 :
142 16 : Simulator::Simulator(MeasurementSet &theMs)
143 16 : : msname_p(""), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
144 16 : ac_p(0), distance_p(0), ve_p(), vc_p(), vp_p(0), gvp_p(0),
145 16 : sim_p(0),
146 : // epJ_p(0),
147 16 : epJTableName_p(),
148 32 : prtlev_(0)
149 : {
150 32 : 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 0 : : msname_p(""), ms_p(0), vs_p(0), seed_p(11111),
188 0 : ac_p(0), distance_p(0),ve_p(), vc_p(), vp_p(0), gvp_p(0),
189 0 : sim_p(0),
190 : // epJ_p(0),
191 0 : 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 47 : Simulator::~Simulator()
232 : {
233 47 : if (ms_p) {
234 47 : ms_p->relinquishAutoLocks();
235 47 : ms_p->unlock();
236 47 : delete ms_p;
237 : }
238 47 : ms_p = 0;
239 47 : if (mssel_p) {
240 43 : mssel_p->relinquishAutoLocks();
241 43 : mssel_p->unlock();
242 43 : delete mssel_p;
243 : }
244 47 : mssel_p = 0;
245 47 : if (vs_p) {
246 43 : delete vs_p;
247 : }
248 47 : vs_p = 0;
249 :
250 : // Delete all vis-plane calibration corruption terms
251 47 : resetviscal();
252 :
253 : // Delete all im-plane calibration corruption terms
254 47 : resetimcal();
255 :
256 47 : if(sim_p) delete sim_p; sim_p = 0;
257 :
258 47 : if(sm_p) delete sm_p; sm_p = 0;
259 47 : if(ft_p) delete ft_p; ft_p = 0;
260 47 : if(cft_p) delete cft_p; cft_p = 0;
261 :
262 47 : }
263 :
264 :
265 47 : void Simulator::defaults()
266 : {
267 47 : UnitMap::putUser("Pixel", UnitVal(1.0), "Pixel solid angle");
268 47 : UnitMap::putUser("Beam", UnitVal(1.0), "Beam solid angle");
269 47 : gridfunction_p="SF";
270 : // Use half the machine memory as cache. The user can override
271 : // this via the setoptions function().
272 47 : cache_p=(HostInfo::memoryTotal()/8)*1024;
273 :
274 47 : tile_p=16;
275 47 : ftmachine_p="gridft";
276 47 : padding_p=1.3;
277 47 : facets_p=1;
278 47 : sm_p = 0;
279 47 : ft_p = 0;
280 47 : cft_p = 0;
281 47 : vp_p = 0;
282 47 : gvp_p = 0;
283 47 : sim_p = 0;
284 47 : images_p = 0;
285 47 : nmodels_p = 1;
286 : // info for configurations
287 47 : areStationCoordsSet_p = false;
288 47 : telescope_p = "UNSET";
289 47 : nmodels_p = 0;
290 :
291 : // info for fields and schedule:
292 47 : nField=0;
293 47 : sourceName_p.resize(1);
294 47 : sourceName_p[0]="UNSET";
295 47 : calCode_p.resize(1);
296 47 : calCode_p[0]="";
297 47 : sourceDirection_p.resize(1);
298 :
299 : // info for spectral windows
300 47 : nSpw=0;
301 47 : spWindowName_p.resize(1);
302 47 : nChan_p.resize(1);
303 47 : startFreq_p.resize(1);
304 47 : freqInc_p.resize(1);
305 47 : freqRes_p.resize(1);
306 47 : stokesString_p.resize(1);
307 47 : spWindowName_p[0]="UNSET";
308 47 : nChan_p[0]=1;
309 47 : startFreq_p[0]=Quantity(50., "GHz");
310 47 : freqInc_p[0]=Quantity(0.1, "MHz");
311 47 : freqRes_p[0]=Quantity(0.1, "MHz");
312 47 : stokesString_p[0]="RR RL LR LL";
313 :
314 : // feeds
315 47 : feedMode_p = "perfect R L";
316 47 : nFeeds_p = 1;
317 47 : feedsHaveBeenSet = false;
318 47 : feedsInitialized = false;
319 :
320 : // times
321 47 : integrationTime_p = Quantity(10.0, "s");
322 47 : useHourAngle_p=true;
323 47 : refTime_p = MEpoch(Quantity(0.0, "s"), MEpoch::UTC);
324 47 : timesHaveBeenSet_p=false;
325 :
326 : // VP stuff
327 47 : doVP_p=false;
328 47 : doDefaultVP_p = true;
329 :
330 47 : };
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 0 : }
355 :
356 47 : Bool Simulator::resetviscal() {
357 94 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
358 : try {
359 :
360 47 : os << "Resetting all visibility corruption components" << LogIO::POST;
361 :
362 : // The noise term (for now)
363 47 : if(ac_p) delete ac_p; ac_p=0;
364 :
365 : // Delete all VisCals
366 79 : for (uInt i=0;i<vc_p.nelements();++i)
367 32 : if (vc_p[i]) delete vc_p[i];
368 47 : vc_p.resize(0,true);
369 :
370 : // reset the VisEquation (by sending an empty vc_p)
371 47 : 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 0 : }
377 47 : return true;
378 47 : }
379 :
380 :
381 47 : Bool Simulator::resetimcal() {
382 94 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
383 : try {
384 :
385 47 : os << "Reset all image-plane corruption components" << LogIO::POST;
386 :
387 47 : if(vp_p) delete vp_p; vp_p=0;
388 47 : 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 0 : }
397 47 : return true;
398 47 : }
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 0 : }
415 0 : return true;
416 0 : }
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 0 : }
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 0 : }
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 29 : Bool Simulator::settimes(const Quantity& integrationTime,
658 : const Bool useHourAngle,
659 : const MEpoch& refTime)
660 : {
661 :
662 58 : LogIO os(LogOrigin("simulator", "settimes()"));
663 : // Negative integration time crashes casa
664 29 : AlwaysAssert(integrationTime.getValue("s")>=0, AipsError);
665 : try {
666 :
667 29 : integrationTime_p=integrationTime;
668 29 : useHourAngle_p=useHourAngle;
669 29 : refTime_p=refTime;
670 :
671 : os << "Times " << endl
672 29 : << " Integration time " << integrationTime.getValue("s") << "s" << LogIO::POST;
673 29 : if(useHourAngle) {
674 29 : os << " Times will be interpreted as hour angles for first source" << LogIO::POST;
675 : }
676 :
677 29 : sim_p->settimes(integrationTime, useHourAngle, refTime);
678 :
679 29 : timesHaveBeenSet_p=true;
680 :
681 29 : return true;
682 0 : } catch (AipsError x) {
683 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
684 0 : return false;
685 0 : }
686 : return true;
687 :
688 29 : }
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 30 : 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 60 : LogIO os(LogOrigin("Simulator", "setconfig()"));
712 :
713 30 : telescope_p = telname;
714 30 : x_p.resize(x.nelements());
715 30 : x_p = x;
716 30 : y_p.resize(y.nelements());
717 30 : y_p = y;
718 30 : z_p.resize(z.nelements());
719 30 : z_p = z;
720 30 : diam_p.resize(dishDiameter.nelements());
721 30 : diam_p = dishDiameter;
722 30 : offset_p.resize(offset.nelements());
723 30 : offset_p = offset;
724 30 : mount_p.resize(mount.nelements());
725 30 : mount_p = mount;
726 30 : antName_p.resize(antName.nelements());
727 30 : antName_p = antName;
728 30 : padName_p.resize(padName.nelements());
729 30 : padName_p = padName;
730 30 : coordsystem_p = coordsystem;
731 30 : mRefLocation_p = mRefLocation;
732 :
733 30 : uInt nn = x_p.nelements();
734 :
735 30 : if (diam_p.nelements() == 1) {
736 10 : diam_p.resize(nn);
737 10 : diam_p.set(dishDiameter(0));
738 : }
739 30 : if (mount_p.nelements() == 1) {
740 28 : mount_p.resize(nn);
741 28 : mount_p.set(mount(0));
742 : }
743 30 : if (mount_p.nelements() == 0) {
744 0 : mount_p.resize(nn);
745 0 : mount_p.set("alt-az");
746 : }
747 30 : if (offset_p.nelements() == 1) {
748 28 : offset_p.resize(nn);
749 28 : offset_p.set(offset(0));
750 : }
751 30 : if (offset_p.nelements() == 0) {
752 0 : offset_p.resize(nn);
753 0 : offset_p.set(0.0);
754 : }
755 30 : if (antName_p.nelements() == 1) {
756 10 : antName_p.resize(nn);
757 10 : antName_p.set(antName(0));
758 : }
759 30 : if (antName_p.nelements() == 0) {
760 0 : antName_p.resize(nn);
761 0 : antName_p.set("UNKNOWN");
762 : }
763 30 : if (padName_p.nelements() == 1) {
764 10 : padName_p.resize(nn);
765 10 : padName_p.set(padName(0));
766 : }
767 30 : if (padName_p.nelements() == 0) {
768 0 : padName_p.resize(nn);
769 0 : padName_p.set("UNKNOWN");
770 : }
771 :
772 30 : AlwaysAssert( (nn == y_p.nelements()) , AipsError);
773 30 : AlwaysAssert( (nn == z_p.nelements()) , AipsError);
774 30 : AlwaysAssert( (nn == diam_p.nelements()) , AipsError);
775 30 : AlwaysAssert( (nn == offset_p.nelements()) , AipsError);
776 30 : AlwaysAssert( (nn == mount_p.nelements()) , AipsError);
777 :
778 30 : areStationCoordsSet_p = true;
779 :
780 30 : sim_p->initAnt(telescope_p, x_p, y_p, z_p, diam_p, offset_p, mount_p, antName_p, padName_p,
781 30 : coordsystem_p, mRefLocation_p);
782 :
783 30 : return true;
784 30 : }
785 :
786 :
787 :
788 16 : Bool Simulator::getconfig() {
789 : // get it from NewMSSimulator
790 16 : 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 16 : }
808 :
809 :
810 :
811 490 : Bool Simulator::setfield(const String& sourceName,
812 : const MDirection& sourceDirection,
813 : const String& calCode,
814 : const Quantity& distance)
815 : {
816 980 : LogIO os(LogOrigin("Simulator", "setfield()"));
817 : try {
818 :
819 490 : if (sourceName == "") {
820 0 : os << LogIO::SEVERE << "must provide a source name" << LogIO::POST;
821 0 : return false;
822 : }
823 :
824 490 : nField++;
825 :
826 490 : if (prtlev()>2) os << "nField = " << nField << LogIO::POST;
827 :
828 490 : distance_p.resize(nField,true);
829 490 : distance_p[nField-1]=distance;
830 490 : sourceName_p.resize(nField,true);
831 490 : sourceName_p[nField-1]=sourceName;
832 490 : sourceDirection_p.resize(nField,true);
833 490 : sourceDirection_p[nField-1]=sourceDirection;
834 490 : calCode_p.resize(nField,true);
835 490 : calCode_p[nField-1]=calCode;
836 :
837 490 : 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 0 : }
843 490 : return true;
844 490 : };
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 0 : }
901 : }
902 :
903 :
904 0 : return true;
905 : }
906 :
907 :
908 :
909 30 : 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 60 : LogIO os(LogOrigin("Simulator", "setspwindow()"));
919 : try {
920 30 : if (nChan == 0) {
921 0 : os << LogIO::SEVERE << "must provide nchannels" << LogIO::POST;
922 0 : return false;
923 : }
924 :
925 30 : nSpw++;
926 : #ifdef RI_DEBUG
927 : os << "nspw = " << nSpw << LogIO::POST;
928 : #endif
929 30 : spWindowName_p.resize(nSpw,true);
930 30 : spWindowName_p[nSpw-1] = spwName;
931 30 : nChan_p.resize(nSpw,true);
932 30 : nChan_p[nSpw-1] = nChan;
933 30 : startFreq_p.resize(nSpw,true);
934 30 : startFreq_p[nSpw-1] = freq;
935 30 : freqInc_p.resize(nSpw,true);
936 30 : freqInc_p[nSpw-1] = deltafreq;
937 30 : freqRes_p.resize(nSpw,true);
938 30 : freqRes_p[nSpw-1] = freqresolution;
939 30 : stokesString_p.resize(nSpw,true);
940 30 : 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 30 : sim_p->initSpWindows(spWindowName_p[nSpw-1], nChan_p[nSpw-1],
948 30 : startFreq_p[nSpw-1], freqInc_p[nSpw-1],
949 30 : freqRes_p[nSpw-1], freqType,
950 30 : 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 0 : }
956 30 : return true;
957 30 : };
958 :
959 :
960 30 : Bool Simulator::setfeed(const String& mode,
961 : const Vector<Double>& x,
962 : const Vector<Double>& y,
963 : const Vector<String>& pol)
964 : {
965 60 : LogIO os(LogOrigin("Simulator", "setfeed()"));
966 :
967 30 : 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 30 : feedMode_p = mode;
974 30 : sim_p->initFeeds(feedMode_p, x, y, pol);
975 :
976 30 : nFeeds_p = x.nelements();
977 30 : feedsHaveBeenSet = true;
978 :
979 30 : return true;
980 30 : };
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 17 : return true;
1022 17 : };
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 32 : 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 16 : 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 16 : 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 16 : 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 16 : } catch (AipsError x) {
1265 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1266 0 : return false;
1267 0 : }
1268 16 : return true;
1269 16 : }
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 0 : } 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 0 : }
1334 0 : return true;
1335 0 : }
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 0 : } 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 0 : }
1454 : #endif
1455 0 : return true;
1456 0 : }
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 0 : }
1526 : #endif
1527 0 : return true;
1528 0 : }
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 0 : }
1596 : return true;
1597 :
1598 0 : }
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 0 : }
1629 : return true;
1630 0 : };
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 0 : }
1661 : return true;
1662 0 : }
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 0 : }
1690 :
1691 : applyPointingOffsets_p = applyPointingOffsets;
1692 : doPBCorrection_p = doPBCorrection;
1693 : return true;
1694 0 : }
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 32 : 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 32 : 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 32 : } 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 0 : }
1745 32 : return true;
1746 32 : }
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 0 : }
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 0 : }
1854 : return false;
1855 0 : }
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 16 : 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 16 : 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 16 : 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 16 : } 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 0 : }
1993 16 : return true;
1994 16 : }
1995 :
1996 :
1997 :
1998 :
1999 :
2000 :
2001 :
2002 :
2003 :
2004 :
2005 :
2006 :
2007 :
2008 :
2009 :
2010 :
2011 :
2012 :
2013 :
2014 :
2015 4 : 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 8 : LogIO os(LogOrigin("Simulator", "observe()", WHERE));
2030 :
2031 :
2032 : try {
2033 :
2034 4 : 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 4 : 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 4 : 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 4 : if(ms_p) delete ms_p; ms_p=0;
2050 4 : if(mssel_p) delete mssel_p; mssel_p=0;
2051 8 : ms_p = new MeasurementSet(msname_p,
2052 4 : TableLock(TableLock::AutoNoReadLocking),
2053 4 : Table::Update);
2054 :
2055 4 : ms_p->flush();
2056 4 : 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 0 : }
2062 4 : return true;
2063 4 : }
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 54 : 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 0 : }
2115 27 : return true;
2116 27 : }
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 54 : 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 27 : 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 : 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 27 : } 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 0 : }
2214 27 : return true;
2215 27 : }
2216 :
2217 :
2218 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2219 : // copied from SynthesisImager
2220 27 : void Simulator::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2221 : {
2222 54 : 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 54 : 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 0 : }
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 23 : }
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 27 : MSColumns msc(*ams);
2381 27 : String telescop=msc.observation().telescopeName()(0);
2382 : PBMath::CommonPB kpb;
2383 27 : 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 10 : }
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 0 : }
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 0 : }
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 0 : }
2563 17 : vp_p->summary();
2564 17 : se_p->setSkyJones(*vp_p);
2565 : } else {
2566 10 : vp_p=0;
2567 : }
2568 27 : } 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 0 : }
2574 27 : return true;
2575 27 : };
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 2 : Bool Simulator::setlimits(const Double shadowLimit,
2603 : const Quantity& elevationLimit)
2604 : {
2605 :
2606 4 : LogIO os(LogOrigin("Simulator", "setlimits()", WHERE));
2607 :
2608 : try {
2609 :
2610 2 : sim_p->setFractionBlockageLimit( shadowLimit );
2611 2 : 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 0 : }
2616 2 : return true;
2617 2 : }
2618 :
2619 30 : Bool Simulator::setauto(const Double autocorrwt)
2620 : {
2621 :
2622 60 : LogIO os(LogOrigin("Simulator", "setauto()", WHERE));
2623 :
2624 : try {
2625 :
2626 30 : sim_p->setAutoCorrelationWt(autocorrwt);
2627 :
2628 : } catch (AipsError x) {
2629 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2630 : return false;
2631 : }
2632 30 : return true;
2633 30 : }
2634 :
2635 129 : void Simulator::makeVisSet() {
2636 :
2637 129 : if(vs_p) return;
2638 :
2639 70 : 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 70 : 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 70 : }
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 10 : return true;
2688 10 : }
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 0 : }
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 0 : }
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 86 : 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 43 : MSDataDescColumns dataDescCol(ms_p->dataDescription());
2737 43 : 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 43 : TableExprNode condition;
2764 43 : String colf=MS::columnName(MS::FIELD_ID);
2765 43 : 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 43 : }
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 0 : } 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 43 : } catch (AipsError x) {
2859 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg()
2860 0 : << LogIO::POST;
2861 0 : return false;
2862 0 : }
2863 : return true;
2864 43 : }
2865 :
2866 : } // end namespace casa
|
