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 26 : Simulator::Simulator(String& msname)
117 26 : : msname_p(msname), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
118 26 : ac_p(0), distance_p(0),ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
119 26 : sim_p(0),
120 : // epJ_p(0),
121 26 : epJTableName_p(),
122 52 : prtlev_(0)
123 : {
124 26 : defaults();
125 :
126 26 : if(!sim_p) {
127 26 : sim_p= new NewMSSimulator(msname);
128 : //sim_p->setPrtlev(prtlev());
129 : }
130 :
131 26 : ve_p.setPrtlev(prtlev());
132 :
133 : // Make a MeasurementSet object for the disk-base MeasurementSet that we just
134 : // created
135 26 : ms_p = new MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
136 26 : Table::Update);
137 26 : AlwaysAssert(ms_p, AipsError);
138 :
139 26 : }
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 42 : Simulator::~Simulator()
232 : {
233 42 : if (ms_p) {
234 42 : ms_p->relinquishAutoLocks();
235 42 : ms_p->unlock();
236 42 : delete ms_p;
237 : }
238 42 : ms_p = 0;
239 42 : if (mssel_p) {
240 38 : mssel_p->relinquishAutoLocks();
241 38 : mssel_p->unlock();
242 38 : delete mssel_p;
243 : }
244 42 : mssel_p = 0;
245 42 : if (vs_p) {
246 38 : delete vs_p;
247 : }
248 42 : vs_p = 0;
249 :
250 : // Delete all vis-plane calibration corruption terms
251 42 : resetviscal();
252 :
253 : // Delete all im-plane calibration corruption terms
254 42 : resetimcal();
255 :
256 42 : if(sim_p) delete sim_p; sim_p = 0;
257 :
258 42 : if(sm_p) delete sm_p; sm_p = 0;
259 42 : if(ft_p) delete ft_p; ft_p = 0;
260 42 : if(cft_p) delete cft_p; cft_p = 0;
261 :
262 42 : }
263 :
264 :
265 42 : void Simulator::defaults()
266 : {
267 42 : UnitMap::putUser("Pixel", UnitVal(1.0), "Pixel solid angle");
268 42 : UnitMap::putUser("Beam", UnitVal(1.0), "Beam solid angle");
269 42 : gridfunction_p="SF";
270 : // Use half the machine memory as cache. The user can override
271 : // this via the setoptions function().
272 42 : cache_p=(HostInfo::memoryTotal()/8)*1024;
273 :
274 42 : tile_p=16;
275 42 : ftmachine_p="gridft";
276 42 : padding_p=1.3;
277 42 : facets_p=1;
278 42 : sm_p = 0;
279 42 : ft_p = 0;
280 42 : cft_p = 0;
281 42 : vp_p = 0;
282 42 : gvp_p = 0;
283 42 : sim_p = 0;
284 42 : images_p = 0;
285 42 : nmodels_p = 1;
286 : // info for configurations
287 42 : areStationCoordsSet_p = false;
288 42 : telescope_p = "UNSET";
289 42 : nmodels_p = 0;
290 :
291 : // info for fields and schedule:
292 42 : nField=0;
293 42 : sourceName_p.resize(1);
294 42 : sourceName_p[0]="UNSET";
295 42 : calCode_p.resize(1);
296 42 : calCode_p[0]="";
297 42 : sourceDirection_p.resize(1);
298 :
299 : // info for spectral windows
300 42 : nSpw=0;
301 42 : spWindowName_p.resize(1);
302 42 : nChan_p.resize(1);
303 42 : startFreq_p.resize(1);
304 42 : freqInc_p.resize(1);
305 42 : freqRes_p.resize(1);
306 42 : stokesString_p.resize(1);
307 42 : spWindowName_p[0]="UNSET";
308 42 : nChan_p[0]=1;
309 42 : startFreq_p[0]=Quantity(50., "GHz");
310 42 : freqInc_p[0]=Quantity(0.1, "MHz");
311 42 : freqRes_p[0]=Quantity(0.1, "MHz");
312 42 : stokesString_p[0]="RR RL LR LL";
313 :
314 : // feeds
315 42 : feedMode_p = "perfect R L";
316 42 : nFeeds_p = 1;
317 42 : feedsHaveBeenSet = false;
318 42 : feedsInitialized = false;
319 :
320 : // times
321 42 : integrationTime_p = Quantity(10.0, "s");
322 42 : useHourAngle_p=true;
323 42 : refTime_p = MEpoch(Quantity(0.0, "s"), MEpoch::UTC);
324 42 : timesHaveBeenSet_p=false;
325 :
326 : // VP stuff
327 42 : doVP_p=false;
328 42 : doDefaultVP_p = true;
329 :
330 42 : };
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 42 : Bool Simulator::resetviscal() {
357 84 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
358 : try {
359 :
360 42 : os << "Resetting all visibility corruption components" << LogIO::POST;
361 :
362 : // The noise term (for now)
363 42 : if(ac_p) delete ac_p; ac_p=0;
364 :
365 : // Delete all VisCals
366 74 : for (uInt i=0;i<vc_p.nelements();++i)
367 32 : if (vc_p[i]) delete vc_p[i];
368 42 : vc_p.resize(0,true);
369 :
370 : // reset the VisEquation (by sending an empty vc_p)
371 42 : 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 42 : return true;
378 42 : }
379 :
380 :
381 42 : Bool Simulator::resetimcal() {
382 84 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
383 : try {
384 :
385 42 : os << "Reset all image-plane corruption components" << LogIO::POST;
386 :
387 42 : if(vp_p) delete vp_p; vp_p=0;
388 42 : 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 42 : return true;
398 42 : }
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 24 : Bool Simulator::settimes(const Quantity& integrationTime,
658 : const Bool useHourAngle,
659 : const MEpoch& refTime)
660 : {
661 :
662 48 : LogIO os(LogOrigin("simulator", "settimes()"));
663 : // Negative integration time crashes casa
664 24 : AlwaysAssert(integrationTime.getValue("s")>=0, AipsError);
665 : try {
666 :
667 24 : integrationTime_p=integrationTime;
668 24 : useHourAngle_p=useHourAngle;
669 24 : refTime_p=refTime;
670 :
671 : os << "Times " << endl
672 24 : << " Integration time " << integrationTime.getValue("s") << "s" << LogIO::POST;
673 24 : if(useHourAngle) {
674 24 : os << " Times will be interpreted as hour angles for first source" << LogIO::POST;
675 : }
676 :
677 24 : sim_p->settimes(integrationTime, useHourAngle, refTime);
678 :
679 24 : timesHaveBeenSet_p=true;
680 :
681 24 : 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 24 : }
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 25 : 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 50 : LogIO os(LogOrigin("Simulator", "setconfig()"));
712 :
713 25 : telescope_p = telname;
714 25 : x_p.resize(x.nelements());
715 25 : x_p = x;
716 25 : y_p.resize(y.nelements());
717 25 : y_p = y;
718 25 : z_p.resize(z.nelements());
719 25 : z_p = z;
720 25 : diam_p.resize(dishDiameter.nelements());
721 25 : diam_p = dishDiameter;
722 25 : offset_p.resize(offset.nelements());
723 25 : offset_p = offset;
724 25 : mount_p.resize(mount.nelements());
725 25 : mount_p = mount;
726 25 : antName_p.resize(antName.nelements());
727 25 : antName_p = antName;
728 25 : padName_p.resize(padName.nelements());
729 25 : padName_p = padName;
730 25 : coordsystem_p = coordsystem;
731 25 : mRefLocation_p = mRefLocation;
732 :
733 25 : uInt nn = x_p.nelements();
734 :
735 25 : if (diam_p.nelements() == 1) {
736 8 : diam_p.resize(nn);
737 8 : diam_p.set(dishDiameter(0));
738 : }
739 25 : if (mount_p.nelements() == 1) {
740 23 : mount_p.resize(nn);
741 23 : mount_p.set(mount(0));
742 : }
743 25 : if (mount_p.nelements() == 0) {
744 0 : mount_p.resize(nn);
745 0 : mount_p.set("alt-az");
746 : }
747 25 : if (offset_p.nelements() == 1) {
748 23 : offset_p.resize(nn);
749 23 : offset_p.set(offset(0));
750 : }
751 25 : if (offset_p.nelements() == 0) {
752 0 : offset_p.resize(nn);
753 0 : offset_p.set(0.0);
754 : }
755 25 : if (antName_p.nelements() == 1) {
756 8 : antName_p.resize(nn);
757 8 : antName_p.set(antName(0));
758 : }
759 25 : if (antName_p.nelements() == 0) {
760 0 : antName_p.resize(nn);
761 0 : antName_p.set("UNKNOWN");
762 : }
763 25 : if (padName_p.nelements() == 1) {
764 8 : padName_p.resize(nn);
765 8 : padName_p.set(padName(0));
766 : }
767 25 : if (padName_p.nelements() == 0) {
768 0 : padName_p.resize(nn);
769 0 : padName_p.set("UNKNOWN");
770 : }
771 :
772 25 : AlwaysAssert( (nn == y_p.nelements()) , AipsError);
773 25 : AlwaysAssert( (nn == z_p.nelements()) , AipsError);
774 25 : AlwaysAssert( (nn == diam_p.nelements()) , AipsError);
775 25 : AlwaysAssert( (nn == offset_p.nelements()) , AipsError);
776 25 : AlwaysAssert( (nn == mount_p.nelements()) , AipsError);
777 :
778 25 : areStationCoordsSet_p = true;
779 :
780 25 : sim_p->initAnt(telescope_p, x_p, y_p, z_p, diam_p, offset_p, mount_p, antName_p, padName_p,
781 25 : coordsystem_p, mRefLocation_p);
782 :
783 25 : return true;
784 25 : }
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 359 : Bool Simulator::setfield(const String& sourceName,
812 : const MDirection& sourceDirection,
813 : const String& calCode,
814 : const Quantity& distance)
815 : {
816 718 : LogIO os(LogOrigin("Simulator", "setfield()"));
817 : try {
818 :
819 359 : if (sourceName == "") {
820 0 : os << LogIO::SEVERE << "must provide a source name" << LogIO::POST;
821 0 : return false;
822 : }
823 :
824 359 : nField++;
825 :
826 359 : if (prtlev()>2) os << "nField = " << nField << LogIO::POST;
827 :
828 359 : distance_p.resize(nField,true);
829 359 : distance_p[nField-1]=distance;
830 359 : sourceName_p.resize(nField,true);
831 359 : sourceName_p[nField-1]=sourceName;
832 359 : sourceDirection_p.resize(nField,true);
833 359 : sourceDirection_p[nField-1]=sourceDirection;
834 359 : calCode_p.resize(nField,true);
835 359 : calCode_p[nField-1]=calCode;
836 :
837 359 : 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 359 : return true;
844 359 : };
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 25 : 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 50 : LogIO os(LogOrigin("Simulator", "setspwindow()"));
919 : try {
920 25 : if (nChan == 0) {
921 0 : os << LogIO::SEVERE << "must provide nchannels" << LogIO::POST;
922 0 : return false;
923 : }
924 :
925 25 : nSpw++;
926 : #ifdef RI_DEBUG
927 : os << "nspw = " << nSpw << LogIO::POST;
928 : #endif
929 25 : spWindowName_p.resize(nSpw,true);
930 25 : spWindowName_p[nSpw-1] = spwName;
931 25 : nChan_p.resize(nSpw,true);
932 25 : nChan_p[nSpw-1] = nChan;
933 25 : startFreq_p.resize(nSpw,true);
934 25 : startFreq_p[nSpw-1] = freq;
935 25 : freqInc_p.resize(nSpw,true);
936 25 : freqInc_p[nSpw-1] = deltafreq;
937 25 : freqRes_p.resize(nSpw,true);
938 25 : freqRes_p[nSpw-1] = freqresolution;
939 25 : stokesString_p.resize(nSpw,true);
940 25 : 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 25 : sim_p->initSpWindows(spWindowName_p[nSpw-1], nChan_p[nSpw-1],
948 25 : startFreq_p[nSpw-1], freqInc_p[nSpw-1],
949 25 : freqRes_p[nSpw-1], freqType,
950 25 : 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 25 : return true;
957 25 : };
958 :
959 :
960 25 : Bool Simulator::setfeed(const String& mode,
961 : const Vector<Double>& x,
962 : const Vector<Double>& y,
963 : const Vector<String>& pol)
964 : {
965 50 : LogIO os(LogOrigin("Simulator", "setfeed()"));
966 :
967 25 : 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 25 : feedMode_p = mode;
974 25 : sim_p->initFeeds(feedMode_p, x, y, pol);
975 :
976 25 : nFeeds_p = x.nelements();
977 25 : feedsHaveBeenSet = true;
978 :
979 25 : return true;
980 25 : };
981 :
982 :
983 :
984 :
985 14 : 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 28 : LogIO os(LogOrigin("Simulator", "setvp()"));
994 :
995 14 : os << "Setting voltage pattern parameters" << LogIO::POST;
996 14 : doVP_p=dovp;
997 14 : doDefaultVP_p = doDefaultVPs;
998 14 : vpTableStr_p = vpTable;
999 14 : if (doSquint) {
1000 14 : squintType_p = BeamSquint::GOFIGURE;
1001 : } else {
1002 0 : squintType_p = BeamSquint::NONE;
1003 : }
1004 14 : parAngleInc_p = parAngleInc;
1005 14 : skyPosThreshold_p = skyPosThreshold;
1006 :
1007 14 : if (doDefaultVP_p) {
1008 0 : os << "Using system default voltage patterns for each telescope" << LogIO::POST;
1009 : } else {
1010 14 : if (vpTableStr_p != String("")) {
1011 0 : os << "Using user defined voltage patterns in Table "<< vpTableStr_p
1012 0 : << LogIO::POST;
1013 : }
1014 : }
1015 14 : if (doSquint) {
1016 14 : os << "Beam Squint will be included in the VP model" << LogIO::POST;
1017 : os << "and the parallactic angle increment is "
1018 14 : << parAngleInc_p.getValue("deg") << " degrees" << LogIO::POST;
1019 : }
1020 14 : pbLimit_p = pbLimit;
1021 14 : return true;
1022 14 : };
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 22 : 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 44 : LogIO os(LogOrigin("Simulator", "observemany()", WHERE));
2084 :
2085 :
2086 : try {
2087 :
2088 22 : 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 22 : 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 22 : 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 22 : if(ms_p) delete ms_p; ms_p=0;
2103 22 : if(mssel_p) delete mssel_p; mssel_p=0;
2104 44 : ms_p = new MeasurementSet(msname_p,
2105 22 : TableLock(TableLock::AutoNoReadLocking),
2106 22 : Table::Update);
2107 :
2108 22 : ms_p->flush();
2109 22 : 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 22 : return true;
2116 22 : }
2117 :
2118 :
2119 :
2120 :
2121 :
2122 :
2123 22 : Bool Simulator::predict(const Vector<String>& modelImage,
2124 : const String& compList,
2125 : const Bool incremental) {
2126 :
2127 44 : 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 22 : " and componentList: " << compList << LogIO::POST;
2137 22 : if (incremental) {
2138 0 : os << "The data column will be incremented" << LogIO::POST;
2139 : } else {
2140 22 : os << "The data column will be replaced" << LogIO::POST;
2141 : }
2142 22 : if(!ms_p) {
2143 : os << "MeasurementSet pointer is null : logic problem!"
2144 0 : << LogIO::EXCEPTION;
2145 : }
2146 :
2147 22 : bool heterogeneous=False;
2148 661 : for (uInt i=1;i<diam_p.nelements();i++)
2149 639 : if (diam_p(i)!=diam_p(0)) heterogeneous=True;
2150 22 : 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 22 : ms_p->lock();
2160 22 : if(mssel_p) mssel_p->lock();
2161 22 : if (!createSkyEquation( modelImage, compList)) {
2162 0 : os << LogIO::SEVERE << "Failed to create SkyEquation" << LogIO::POST;
2163 0 : return false;
2164 : }
2165 22 : if (incremental) {
2166 0 : se_p->predict(false,MS::MODEL_DATA);
2167 : } else {
2168 22 : se_p->predict(false,MS::DATA); //20091030 RI changed SE::predict to use DATA
2169 : }
2170 22 : destroySkyEquation();
2171 :
2172 : // Copy the predicted visibilities over to the observed and
2173 : // the corrected data columns
2174 22 : makeVisSet();
2175 :
2176 22 : VisIter& vi = vs_p->iter();
2177 22 : VisBuffer vb(vi);
2178 22 : vi.origin();
2179 22 : vi.originChunks();
2180 :
2181 : //os << "Copying predicted visibilities from MODEL_DATA to DATA and CORRECTED_DATA" << LogIO::POST;
2182 :
2183 378 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()){
2184 935 : for (vi.origin(); vi.more(); vi++) {
2185 : // vb.setVisCube(vb.modelVisCube());
2186 :
2187 579 : 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 579 : vi.setVis(vb.visCube(),VisibilityIterator::Corrected);
2200 : }
2201 579 : vb.setModelVisCube(Complex(1.0,0.0));
2202 579 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
2203 : }
2204 : }
2205 22 : ms_p->unlock();
2206 22 : if(mssel_p) mssel_p->lock();
2207 :
2208 22 : } 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 22 : return true;
2215 22 : }
2216 :
2217 :
2218 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2219 : // copied from SynthesisImager
2220 22 : void Simulator::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2221 : {
2222 44 : LogIO os(LogOrigin("Simulator", "getVPRecord",WHERE));
2223 :
2224 22 : VPManager *vpman=VPManager::Instance();
2225 22 : if ((doDefaultVP_p)||(vpTableStr_p != String(""))) {
2226 8 : if (doDefaultVP_p) {
2227 8 : os << "Using default Voltage Patterns from the VPManager" << LogIO::POST;
2228 8 : vpman->reset();
2229 : } else {
2230 0 : os << "Loading Voltage Pattern information from " << vpTableStr_p << LogIO::POST;
2231 0 : vpman->loadfromtable(vpTableStr_p );
2232 : }
2233 8 : 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 14 : os << "Using Voltage Patterns from the VPManager" << LogIO::POST;
2239 : }
2240 :
2241 : // PBMath::CommonPB kpb;
2242 22 : PBMath::enumerateCommonPB(telescop, kpb);
2243 : // Record rec;
2244 22 : vpman->getvp(rec, telescop);
2245 :
2246 : //ostringstream oss; rec.print(oss);
2247 : //os << "Using VP model : " << oss.str() << LogIO::POST;
2248 :
2249 22 : 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 22 : }// get VPRecord
2263 :
2264 :
2265 :
2266 22 : Bool Simulator::createSkyEquation(const Vector<String>& image,
2267 : const String complist)
2268 : {
2269 :
2270 44 : LogIO os(LogOrigin("Simulator", "createSkyEquation()", WHERE));
2271 :
2272 : try {
2273 22 : if(sm_p==0) {
2274 22 : sm_p = new CleanImageSkyModel();
2275 : }
2276 22 : AlwaysAssert(sm_p, AipsError);
2277 :
2278 : // Add the componentlist
2279 22 : 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 15 : componentList_p=0;
2298 : }
2299 :
2300 22 : nmodels_p = image.nelements();
2301 22 : if (nmodels_p == 1 && image(0) == "") nmodels_p = 0;
2302 :
2303 22 : int models_found=0;
2304 22 : if (nmodels_p > 0) {
2305 18 : images_p.resize(nmodels_p);
2306 :
2307 36 : for (Int model=0;model<Int(nmodels_p);model++) {
2308 18 : 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 18 : if(!Table::isReadable(image(model))) {
2314 0 : os << LogIO::SEVERE << image(model) << " is unreadable" << LogIO::POST;
2315 : } else {
2316 18 : images_p[model]=0;
2317 : os << "Opening model " << model << " named "
2318 18 : << image(model) << LogIO::POST;
2319 18 : images_p[model]=new PagedImage<Float>(image(model));
2320 18 : AlwaysAssert(images_p[model], AipsError);
2321 : // RI TODO is this a logic problem with more than one source??
2322 : // Add distance
2323 18 : 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 18 : CoordinateSystem cs = images_p[model]->coordinates();
2333 18 : String errorMsg;
2334 18 : cs.setSpectralConversion(errorMsg,MFrequency::showType(MFrequency::LSRK));
2335 18 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
2336 18 : AlwaysAssert(spectralIndex>=0, AipsError);
2337 18 : SpectralCoordinate spectralCoord=cs.spectralCoordinate(spectralIndex);
2338 18 : Vector<String> units(1); units = "Hz";
2339 : // doesn't work: cs.spectralCoordinate(spectralIndex).setWorldAxisUnits(units);
2340 18 : spectralCoord.setWorldAxisUnits(units);
2341 : // put spectralCoord back into cs
2342 18 : cs.replaceCoordinate(spectralCoord,spectralIndex);
2343 : // and cs into image
2344 18 : images_p[model]->setCoordinateInfo(cs);
2345 :
2346 :
2347 18 : 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 18 : models_found++;
2352 18 : }
2353 : }
2354 : }
2355 : }
2356 :
2357 22 : 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 22 : if(vs_p) {
2365 22 : delete vs_p; vs_p=0;
2366 : }
2367 22 : makeVisSet();
2368 :
2369 22 : cft_p = new SimpleComponentFTMachine();
2370 :
2371 22 : MeasurementSet *ams=0;
2372 :
2373 22 : if(mssel_p) {
2374 22 : ams=mssel_p;
2375 : }
2376 : else {
2377 0 : ams=ms_p;
2378 : }
2379 : // 2016 from SynthesisImager:
2380 22 : MSColumns msc(*ams);
2381 22 : String telescop=msc.observation().telescopeName()(0);
2382 : PBMath::CommonPB kpb;
2383 22 : Record rec;
2384 22 : getVPRecord( rec, kpb, telescop );
2385 :
2386 22 : if((ftmachine_p=="sd")||(ftmachine_p=="both")||(ftmachine_p=="mosaic")) {
2387 8 : if(!gvp_p) {
2388 8 : if (rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2389 8 : String commonPBName;
2390 8 : PBMath::nameCommonPB(kpb, commonPBName);
2391 8 : os << "Using common PB "<<commonPBName<<" for beam calculation for telescope " << telescop << LogIO::POST;
2392 8 : gvp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p, skyPosThreshold_p);
2393 8 : }
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 8 : if(ftmachine_p=="sd") {
2404 8 : os << "Performing Single dish gridding " << LogIO::POST;
2405 8 : if(gridfunction_p=="pb") {
2406 8 : 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 8 : 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 8 : vi.setRowBlocking(100);
2454 : }
2455 :
2456 : else {
2457 14 : os << "Synthesis gridding " << LogIO::POST;
2458 : // Now make the FTMachine
2459 : // if(wprojPlanes_p>1) {
2460 14 : 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 14 : 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 14 : os << "Fourier transforms will use image centers as tangent points" << LogIO::POST;
2535 14 : ft_p = new GridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p, padding_p);
2536 : }
2537 : }
2538 22 : 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 22 : 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 22 : if(doVP_p) {
2550 14 : 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 14 : 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 14 : vp_p->summary();
2564 14 : se_p->setSkyJones(*vp_p);
2565 : } else {
2566 8 : vp_p=0;
2567 : }
2568 22 : } 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 22 : return true;
2575 22 : };
2576 :
2577 22 : void Simulator::destroySkyEquation()
2578 : {
2579 22 : if(se_p) delete se_p; se_p=0;
2580 22 : if(sm_p) delete sm_p; sm_p=0;
2581 22 : if(vp_p) delete vp_p; vp_p=0;
2582 22 : if(componentList_p) delete componentList_p; componentList_p=0;
2583 :
2584 40 : for (Int model=0;model<Int(nmodels_p);model++) {
2585 18 : if(images_p[model]) delete images_p[model]; images_p[model]=0;
2586 : }
2587 22 : };
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 25 : Bool Simulator::setauto(const Double autocorrwt)
2620 : {
2621 :
2622 50 : LogIO os(LogOrigin("Simulator", "setauto()", WHERE));
2623 :
2624 : try {
2625 :
2626 25 : sim_p->setAutoCorrelationWt(autocorrwt);
2627 :
2628 : } catch (AipsError x) {
2629 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2630 : return false;
2631 : }
2632 25 : return true;
2633 25 : }
2634 :
2635 114 : void Simulator::makeVisSet() {
2636 :
2637 114 : if(vs_p) return;
2638 :
2639 60 : Block<int> sort(0);
2640 60 : sort.resize(5);
2641 60 : sort[0] = MS::FIELD_ID;
2642 60 : sort[1] = MS::FEED1;
2643 60 : sort[2] = MS::ARRAY_ID;
2644 60 : sort[3] = MS::DATA_DESC_ID;
2645 60 : sort[4] = MS::TIME;
2646 60 : Matrix<Int> noselection;
2647 60 : if(mssel_p) {
2648 60 : 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 60 : AlwaysAssert(vs_p, AipsError);
2657 :
2658 60 : }
2659 :
2660 :
2661 :
2662 8 : 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 16 : LogIO os(LogOrigin("Simulator", "setoptions()", WHERE));
2668 :
2669 8 : os << "Setting processing options" << LogIO::POST;
2670 :
2671 8 : sim_p->setMaxData(maxData*1024.0*1024.0);
2672 :
2673 8 : ftmachine_p=downcase(ftmachine);
2674 8 : if(cache>0) cache_p=cache;
2675 8 : if(tile>0) tile_p=tile;
2676 8 : gridfunction_p=downcase(gridfunction);
2677 8 : mLocation_p=mLocation;
2678 8 : if(padding>=1.0) {
2679 8 : padding_p=padding;
2680 : }
2681 8 : facets_p=facets;
2682 8 : wprojPlanes_p = wprojPlanes;
2683 : // Destroy the FTMachine
2684 8 : if(ft_p) {delete ft_p; ft_p=0;}
2685 8 : if(cft_p) {delete cft_p; cft_p=0;}
2686 :
2687 8 : return true;
2688 8 : }
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 38 : Bool Simulator::setdata(const Vector<Int>& spectralwindowids,
2717 : const Vector<Int>& fieldids,
2718 : const String& msSelect)
2719 :
2720 : {
2721 :
2722 :
2723 76 : LogIO os(LogOrigin("Simulator", "setdata()", WHERE));
2724 :
2725 38 : 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 38 : os << "Selecting data" << LogIO::POST;
2734 :
2735 : // Map the selected spectral window ids to data description ids
2736 38 : MSDataDescColumns dataDescCol(ms_p->dataDescription());
2737 38 : Vector<Int> ddSpwIds=dataDescCol.spectralWindowId().getColumn();
2738 :
2739 38 : Vector<Int> datadescids(0);
2740 76 : for (uInt row=0; row<ddSpwIds.nelements(); row++) {
2741 38 : Bool found=false;
2742 76 : for (uInt j=0; j<spectralwindowids.nelements(); j++) {
2743 38 : if (ddSpwIds(row)==spectralwindowids(j)) found=true;
2744 : };
2745 38 : if (found) {
2746 38 : datadescids.resize(datadescids.nelements()+1,true);
2747 38 : datadescids(datadescids.nelements()-1)=row;
2748 : };
2749 : };
2750 :
2751 38 : if(vs_p) delete vs_p; vs_p=0;
2752 38 : 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 38 : if(fieldids.nelements()>0||datadescids.nelements()>0) {
2758 38 : os << "Performing selection on MeasurementSet" << LogIO::POST;
2759 38 : Table& original=*ms_p;
2760 :
2761 : // Now we make a condition to do the old FIELD_ID, SPECTRAL_WINDOW_ID
2762 : // selection
2763 38 : TableExprNode condition;
2764 38 : String colf=MS::columnName(MS::FIELD_ID);
2765 38 : String cols=MS::columnName(MS::DATA_DESC_ID);
2766 38 : if(fieldids.nelements()>0&&datadescids.nelements()>0){
2767 22 : condition=original.col(colf).in(fieldids)&&original.col(cols).in(datadescids);
2768 22 : 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 38 : mssel_p = new MeasurementSet(original(condition));
2781 :
2782 : //AlwaysAssert(mssel_p, AipsError);
2783 : //mssel_p->rename(msname_p+"/SELECTED_TABLE", Table::Scratch);
2784 38 : 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 38 : mssel_p->flush();
2793 : }
2794 :
2795 38 : }
2796 : else {
2797 0 : mssel_p=new MeasurementSet(*ms_p);
2798 : }
2799 : {
2800 38 : Int len = msSelect.length();
2801 38 : Int nspace = msSelect.freq (' ');
2802 38 : Bool nullSelect=(msSelect.empty() || nspace==len);
2803 38 : 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 38 : os << "No selection string given" << LogIO::POST;
2827 : }
2828 :
2829 38 : 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 37 : os << "Selection did not drop any rows" << LogIO::POST;
2835 : }
2836 : }
2837 :
2838 : // Now create the VisSet
2839 38 : if(vs_p) delete vs_p; vs_p=0;
2840 38 : makeVisSet();
2841 : //Now assign the source directions to something selected or sensible
2842 : {
2843 : //Int fieldsel=0;
2844 38 : if(fieldids.nelements() >0) {
2845 : //fieldsel=fieldids(0);
2846 : // RI TODO does sim:setdata need this?
2847 22 : nField=fieldids.nelements();
2848 378 : 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 356 : (vs_p->iter()).msColumns().field().name().get(i,sourceName_p[i]);
2852 356 : sourceDirection_p[i]=(vs_p->iter()).msColumns().field().phaseDirMeas(i);
2853 356 : (vs_p->iter()).msColumns().field().code().get(i,calCode_p[i]);
2854 : }
2855 : }
2856 : }
2857 38 : return true;
2858 38 : } 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 38 : }
2865 :
2866 : } // end namespace casa
|
