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 2 : Simulator::Simulator(String& msname)
117 2 : : msname_p(msname), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
118 2 : ac_p(0), distance_p(0),ve_p(), vc_p(), nSpw(0), vp_p(0), gvp_p(0),
119 2 : sim_p(0),
120 : // epJ_p(0),
121 2 : epJTableName_p(),
122 4 : prtlev_(0)
123 : {
124 2 : defaults();
125 :
126 2 : if(!sim_p) {
127 2 : sim_p= new NewMSSimulator(msname);
128 : //sim_p->setPrtlev(prtlev());
129 : }
130 :
131 2 : ve_p.setPrtlev(prtlev());
132 :
133 : // Make a MeasurementSet object for the disk-base MeasurementSet that we just
134 : // created
135 2 : ms_p = new MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
136 2 : Table::Update);
137 2 : AlwaysAssert(ms_p, AipsError);
138 :
139 2 : }
140 :
141 :
142 0 : Simulator::Simulator(MeasurementSet &theMs)
143 0 : : msname_p(""), ms_p(0), mssel_p(0), vs_p(0), seed_p(11111),
144 0 : ac_p(0), distance_p(0), ve_p(), vc_p(), vp_p(0), gvp_p(0),
145 0 : sim_p(0),
146 : // epJ_p(0),
147 0 : epJTableName_p(),
148 0 : prtlev_(0)
149 : {
150 0 : LogIO os(LogOrigin("simulator", "simulator(MeasurementSet& theMs)"));
151 :
152 0 : defaults();
153 :
154 0 : msname_p=theMs.tableName();
155 :
156 0 : if(!sim_p) {
157 0 : sim_p= new NewMSSimulator(theMs);
158 : //sim_p->setPrtlev(prtlev());
159 : }
160 :
161 0 : ve_p.setPrtlev(prtlev());
162 :
163 0 : ms_p = new MeasurementSet(theMs);
164 0 : AlwaysAssert(ms_p, AipsError);
165 :
166 : // get info from the MS into Simulator:
167 0 : if (!getconfig())
168 0 : os << "Can't find antenna information for loaded MS" << LogIO::WARN;
169 0 : 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 0 : freqRes_p.resize(nSpw);
172 0 : for (Int i=0;i<nSpw;++i)
173 0 : freqRes_p[i]=freqInc_p[i];
174 0 : 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 0 : if (!sim_p->getFeedMode(feedMode_p))
178 0 : os << "Can't find Feed information for loaded MS" << LogIO::WARN;
179 : else
180 0 : feedsHaveBeenSet=true;
181 :
182 0 : }
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 2 : Simulator::~Simulator()
232 : {
233 2 : if (ms_p) {
234 2 : ms_p->relinquishAutoLocks();
235 2 : ms_p->unlock();
236 2 : delete ms_p;
237 : }
238 2 : ms_p = 0;
239 2 : if (mssel_p) {
240 0 : mssel_p->relinquishAutoLocks();
241 0 : mssel_p->unlock();
242 0 : delete mssel_p;
243 : }
244 2 : mssel_p = 0;
245 2 : if (vs_p) {
246 0 : delete vs_p;
247 : }
248 2 : vs_p = 0;
249 :
250 : // Delete all vis-plane calibration corruption terms
251 2 : resetviscal();
252 :
253 : // Delete all im-plane calibration corruption terms
254 2 : resetimcal();
255 :
256 2 : if(sim_p) delete sim_p; sim_p = 0;
257 :
258 2 : if(sm_p) delete sm_p; sm_p = 0;
259 2 : if(ft_p) delete ft_p; ft_p = 0;
260 2 : if(cft_p) delete cft_p; cft_p = 0;
261 :
262 2 : }
263 :
264 :
265 2 : void Simulator::defaults()
266 : {
267 2 : UnitMap::putUser("Pixel", UnitVal(1.0), "Pixel solid angle");
268 2 : UnitMap::putUser("Beam", UnitVal(1.0), "Beam solid angle");
269 2 : gridfunction_p="SF";
270 : // Use half the machine memory as cache. The user can override
271 : // this via the setoptions function().
272 2 : cache_p=(HostInfo::memoryTotal()/8)*1024;
273 :
274 2 : tile_p=16;
275 2 : ftmachine_p="gridft";
276 2 : padding_p=1.3;
277 2 : facets_p=1;
278 2 : sm_p = 0;
279 2 : ft_p = 0;
280 2 : cft_p = 0;
281 2 : vp_p = 0;
282 2 : gvp_p = 0;
283 2 : sim_p = 0;
284 2 : images_p = 0;
285 2 : nmodels_p = 1;
286 : // info for configurations
287 2 : areStationCoordsSet_p = false;
288 2 : telescope_p = "UNSET";
289 2 : nmodels_p = 0;
290 :
291 : // info for fields and schedule:
292 2 : nField=0;
293 2 : sourceName_p.resize(1);
294 2 : sourceName_p[0]="UNSET";
295 2 : calCode_p.resize(1);
296 2 : calCode_p[0]="";
297 2 : sourceDirection_p.resize(1);
298 :
299 : // info for spectral windows
300 2 : nSpw=0;
301 2 : spWindowName_p.resize(1);
302 2 : nChan_p.resize(1);
303 2 : startFreq_p.resize(1);
304 2 : freqInc_p.resize(1);
305 2 : freqRes_p.resize(1);
306 2 : stokesString_p.resize(1);
307 2 : spWindowName_p[0]="UNSET";
308 2 : nChan_p[0]=1;
309 2 : startFreq_p[0]=Quantity(50., "GHz");
310 2 : freqInc_p[0]=Quantity(0.1, "MHz");
311 2 : freqRes_p[0]=Quantity(0.1, "MHz");
312 2 : stokesString_p[0]="RR RL LR LL";
313 :
314 : // feeds
315 2 : feedMode_p = "perfect R L";
316 2 : nFeeds_p = 1;
317 2 : feedsHaveBeenSet = false;
318 2 : feedsInitialized = false;
319 :
320 : // times
321 2 : integrationTime_p = Quantity(10.0, "s");
322 2 : useHourAngle_p=true;
323 2 : refTime_p = MEpoch(Quantity(0.0, "s"), MEpoch::UTC);
324 2 : timesHaveBeenSet_p=false;
325 :
326 : // VP stuff
327 2 : doVP_p=false;
328 2 : doDefaultVP_p = true;
329 :
330 2 : };
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 2 : Bool Simulator::resetviscal() {
357 4 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
358 : try {
359 :
360 2 : os << "Resetting all visibility corruption components" << LogIO::POST;
361 :
362 : // The noise term (for now)
363 2 : if(ac_p) delete ac_p; ac_p=0;
364 :
365 : // Delete all VisCals
366 2 : for (uInt i=0;i<vc_p.nelements();++i)
367 0 : if (vc_p[i]) delete vc_p[i];
368 2 : vc_p.resize(0,true);
369 :
370 : // reset the VisEquation (by sending an empty vc_p)
371 2 : 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 2 : return true;
378 2 : }
379 :
380 :
381 2 : Bool Simulator::resetimcal() {
382 4 : LogIO os(LogOrigin("Simulator", "reset()", WHERE));
383 : try {
384 :
385 2 : os << "Reset all image-plane corruption components" << LogIO::POST;
386 :
387 2 : if(vp_p) delete vp_p; vp_p=0;
388 2 : 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 2 : return true;
398 2 : }
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 2 : Bool Simulator::settimes(const Quantity& integrationTime,
658 : const Bool useHourAngle,
659 : const MEpoch& refTime)
660 : {
661 :
662 4 : LogIO os(LogOrigin("simulator", "settimes()"));
663 : // Negative integration time crashes casa
664 2 : AlwaysAssert(integrationTime.getValue("s")>=0, AipsError);
665 : try {
666 :
667 2 : integrationTime_p=integrationTime;
668 2 : useHourAngle_p=useHourAngle;
669 2 : refTime_p=refTime;
670 :
671 : os << "Times " << endl
672 2 : << " Integration time " << integrationTime.getValue("s") << "s" << LogIO::POST;
673 2 : if(useHourAngle) {
674 2 : os << " Times will be interpreted as hour angles for first source" << LogIO::POST;
675 : }
676 :
677 2 : sim_p->settimes(integrationTime, useHourAngle, refTime);
678 :
679 2 : timesHaveBeenSet_p=true;
680 :
681 2 : 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 2 : }
689 :
690 :
691 :
692 0 : Bool Simulator::setseed(const Int seed) {
693 0 : seed_p = seed;
694 0 : return true;
695 : }
696 :
697 :
698 :
699 2 : 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 4 : LogIO os(LogOrigin("Simulator", "setconfig()"));
712 :
713 2 : telescope_p = telname;
714 2 : x_p.resize(x.nelements());
715 2 : x_p = x;
716 2 : y_p.resize(y.nelements());
717 2 : y_p = y;
718 2 : z_p.resize(z.nelements());
719 2 : z_p = z;
720 2 : diam_p.resize(dishDiameter.nelements());
721 2 : diam_p = dishDiameter;
722 2 : offset_p.resize(offset.nelements());
723 2 : offset_p = offset;
724 2 : mount_p.resize(mount.nelements());
725 2 : mount_p = mount;
726 2 : antName_p.resize(antName.nelements());
727 2 : antName_p = antName;
728 2 : padName_p.resize(padName.nelements());
729 2 : padName_p = padName;
730 2 : coordsystem_p = coordsystem;
731 2 : mRefLocation_p = mRefLocation;
732 :
733 2 : uInt nn = x_p.nelements();
734 :
735 2 : if (diam_p.nelements() == 1) {
736 0 : diam_p.resize(nn);
737 0 : diam_p.set(dishDiameter(0));
738 : }
739 2 : if (mount_p.nelements() == 1) {
740 0 : mount_p.resize(nn);
741 0 : mount_p.set(mount(0));
742 : }
743 2 : if (mount_p.nelements() == 0) {
744 0 : mount_p.resize(nn);
745 0 : mount_p.set("alt-az");
746 : }
747 2 : if (offset_p.nelements() == 1) {
748 0 : offset_p.resize(nn);
749 0 : offset_p.set(offset(0));
750 : }
751 2 : if (offset_p.nelements() == 0) {
752 0 : offset_p.resize(nn);
753 0 : offset_p.set(0.0);
754 : }
755 2 : if (antName_p.nelements() == 1) {
756 0 : antName_p.resize(nn);
757 0 : antName_p.set(antName(0));
758 : }
759 2 : if (antName_p.nelements() == 0) {
760 0 : antName_p.resize(nn);
761 0 : antName_p.set("UNKNOWN");
762 : }
763 2 : if (padName_p.nelements() == 1) {
764 0 : padName_p.resize(nn);
765 0 : padName_p.set(padName(0));
766 : }
767 2 : if (padName_p.nelements() == 0) {
768 0 : padName_p.resize(nn);
769 0 : padName_p.set("UNKNOWN");
770 : }
771 :
772 2 : AlwaysAssert( (nn == y_p.nelements()) , AipsError);
773 2 : AlwaysAssert( (nn == z_p.nelements()) , AipsError);
774 2 : AlwaysAssert( (nn == diam_p.nelements()) , AipsError);
775 2 : AlwaysAssert( (nn == offset_p.nelements()) , AipsError);
776 2 : AlwaysAssert( (nn == mount_p.nelements()) , AipsError);
777 :
778 2 : areStationCoordsSet_p = true;
779 :
780 2 : sim_p->initAnt(telescope_p, x_p, y_p, z_p, diam_p, offset_p, mount_p, antName_p, padName_p,
781 2 : coordsystem_p, mRefLocation_p);
782 :
783 2 : return true;
784 2 : }
785 :
786 :
787 :
788 0 : Bool Simulator::getconfig() {
789 : // get it from NewMSSimulator
790 0 : Matrix<Double> xyz_p;
791 : Int nAnt;
792 0 : if (sim_p->getAnt(telescope_p, nAnt, &xyz_p, diam_p, offset_p, mount_p, antName_p, padName_p,
793 0 : coordsystem_p, mRefLocation_p)) {
794 0 : x_p.resize(nAnt);
795 0 : y_p.resize(nAnt);
796 0 : z_p.resize(nAnt);
797 0 : for (Int i=0;i<nAnt;i++) {
798 0 : x_p(i)=xyz_p(0,i);
799 0 : y_p(i)=xyz_p(1,i);
800 0 : z_p(i)=xyz_p(2,i);
801 : }
802 0 : areStationCoordsSet_p = true;
803 0 : return true;
804 : } else {
805 0 : return false;
806 : }
807 0 : }
808 :
809 :
810 :
811 2 : Bool Simulator::setfield(const String& sourceName,
812 : const MDirection& sourceDirection,
813 : const String& calCode,
814 : const Quantity& distance)
815 : {
816 4 : LogIO os(LogOrigin("Simulator", "setfield()"));
817 : try {
818 :
819 2 : if (sourceName == "") {
820 0 : os << LogIO::SEVERE << "must provide a source name" << LogIO::POST;
821 0 : return false;
822 : }
823 :
824 2 : nField++;
825 :
826 2 : if (prtlev()>2) os << "nField = " << nField << LogIO::POST;
827 :
828 2 : distance_p.resize(nField,true);
829 2 : distance_p[nField-1]=distance;
830 2 : sourceName_p.resize(nField,true);
831 2 : sourceName_p[nField-1]=sourceName;
832 2 : sourceDirection_p.resize(nField,true);
833 2 : sourceDirection_p[nField-1]=sourceDirection;
834 2 : calCode_p.resize(nField,true);
835 2 : calCode_p[nField-1]=calCode;
836 :
837 2 : 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 2 : return true;
844 2 : };
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 2 : 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 4 : LogIO os(LogOrigin("Simulator", "setspwindow()"));
919 : try {
920 2 : if (nChan == 0) {
921 0 : os << LogIO::SEVERE << "must provide nchannels" << LogIO::POST;
922 0 : return false;
923 : }
924 :
925 2 : nSpw++;
926 : #ifdef RI_DEBUG
927 : os << "nspw = " << nSpw << LogIO::POST;
928 : #endif
929 2 : spWindowName_p.resize(nSpw,true);
930 2 : spWindowName_p[nSpw-1] = spwName;
931 2 : nChan_p.resize(nSpw,true);
932 2 : nChan_p[nSpw-1] = nChan;
933 2 : startFreq_p.resize(nSpw,true);
934 2 : startFreq_p[nSpw-1] = freq;
935 2 : freqInc_p.resize(nSpw,true);
936 2 : freqInc_p[nSpw-1] = deltafreq;
937 2 : freqRes_p.resize(nSpw,true);
938 2 : freqRes_p[nSpw-1] = freqresolution;
939 2 : stokesString_p.resize(nSpw,true);
940 2 : 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 2 : sim_p->initSpWindows(spWindowName_p[nSpw-1], nChan_p[nSpw-1],
948 2 : startFreq_p[nSpw-1], freqInc_p[nSpw-1],
949 2 : freqRes_p[nSpw-1], freqType,
950 2 : 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 2 : return true;
957 2 : };
958 :
959 :
960 2 : Bool Simulator::setfeed(const String& mode,
961 : const Vector<Double>& x,
962 : const Vector<Double>& y,
963 : const Vector<String>& pol)
964 : {
965 4 : LogIO os(LogOrigin("Simulator", "setfeed()"));
966 :
967 2 : 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 2 : feedMode_p = mode;
974 2 : sim_p->initFeeds(feedMode_p, x, y, pol);
975 :
976 2 : nFeeds_p = x.nelements();
977 2 : feedsHaveBeenSet = true;
978 :
979 2 : return true;
980 2 : };
981 :
982 :
983 :
984 :
985 0 : 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 0 : LogIO os(LogOrigin("Simulator", "setvp()"));
994 :
995 0 : os << "Setting voltage pattern parameters" << LogIO::POST;
996 0 : doVP_p=dovp;
997 0 : doDefaultVP_p = doDefaultVPs;
998 0 : vpTableStr_p = vpTable;
999 0 : if (doSquint) {
1000 0 : squintType_p = BeamSquint::GOFIGURE;
1001 : } else {
1002 0 : squintType_p = BeamSquint::NONE;
1003 : }
1004 0 : parAngleInc_p = parAngleInc;
1005 0 : skyPosThreshold_p = skyPosThreshold;
1006 :
1007 0 : if (doDefaultVP_p) {
1008 0 : os << "Using system default voltage patterns for each telescope" << LogIO::POST;
1009 : } else {
1010 0 : if (vpTableStr_p != String("")) {
1011 0 : os << "Using user defined voltage patterns in Table "<< vpTableStr_p
1012 0 : << LogIO::POST;
1013 : }
1014 : }
1015 0 : if (doSquint) {
1016 0 : os << "Beam Squint will be included in the VP model" << LogIO::POST;
1017 : os << "and the parallactic angle increment is "
1018 0 : << parAngleInc_p.getValue("deg") << " degrees" << LogIO::POST;
1019 : }
1020 0 : pbLimit_p = pbLimit;
1021 0 : return true;
1022 0 : };
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 0 : 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 0 : LogIO os(LogOrigin("Simulator", "setnoise2()", WHERE));
1067 :
1068 : try {
1069 :
1070 : //cout<<" Sim::setnoise "<<pground<<" "<<relhum<<" "<<altitude<<" "<<waterheight<<" "<<pwv<<endl;
1071 :
1072 0 : noisemode_p = mode;
1073 :
1074 0 : RecordDesc simparDesc;
1075 0 : simparDesc.addField ("type", TpString);
1076 0 : simparDesc.addField ("mode", TpString);
1077 0 : simparDesc.addField ("caltable", TpString);
1078 :
1079 0 : simparDesc.addField ("amplitude", TpFloat); // for constant scale
1080 : // simparDesc.addField ("scale", TpFloat); // for fractional fluctuations
1081 0 : simparDesc.addField ("combine" ,TpString);
1082 :
1083 : // have to be defined here or else have to be added later
1084 0 : simparDesc.addField ("startTime", TpDouble);
1085 0 : simparDesc.addField ("stopTime", TpDouble);
1086 :
1087 0 : simparDesc.addField ("antefficiency" ,TpFloat);
1088 0 : simparDesc.addField ("spillefficiency",TpFloat);
1089 0 : simparDesc.addField ("correfficiency" ,TpFloat);
1090 0 : simparDesc.addField ("trx" ,TpFloat);
1091 0 : simparDesc.addField ("tground" ,TpFloat);
1092 0 : simparDesc.addField ("tcmb" ,TpFloat);
1093 0 : simparDesc.addField ("senscoeff" ,TpFloat);
1094 0 : simparDesc.addField ("rxType" ,TpInt);
1095 :
1096 : // user-override of ATM calculated tau
1097 0 : simparDesc.addField ("tatmos" ,TpFloat);
1098 0 : simparDesc.addField ("tau0" ,TpFloat);
1099 :
1100 0 : simparDesc.addField ("mean_pwv" ,TpDouble);
1101 0 : simparDesc.addField ("pground" ,TpDouble);
1102 0 : simparDesc.addField ("relhum" ,TpFloat);
1103 0 : simparDesc.addField ("altitude" ,TpDouble);
1104 0 : simparDesc.addField ("waterheight" ,TpDouble);
1105 :
1106 0 : simparDesc.addField ("seed" ,TpInt);
1107 :
1108 : // RI todo setnoise2 if tau0 is not defined, use freqdep
1109 :
1110 0 : String caltbl=caltable;
1111 0 : caltbl.trim();
1112 : string::size_type strlen;
1113 0 : strlen=caltbl.length();
1114 0 : if (strlen>3)
1115 0 : if (caltbl.substr(strlen-3,3)=="cal") {
1116 0 : caltbl.resize(strlen-3);
1117 0 : strlen-=3;
1118 : }
1119 0 : if (strlen>1)
1120 0 : if (caltbl.substr(strlen-1,1)==".") {
1121 0 : caltbl.resize(strlen-1);
1122 0 : strlen-=1;
1123 : }
1124 0 : if (strlen>1)
1125 0 : if (caltbl.substr(strlen-1,1)=="_") {
1126 0 : caltbl.resize(strlen-1);
1127 0 : strlen-=1;
1128 : }
1129 :
1130 0 : Record simpar(simparDesc,RecordInterface::Variable);
1131 0 : simpar.define ("seed", seed_p);
1132 0 : simpar.define ("type", "A Noise");
1133 0 : if (strlen>1)
1134 0 : simpar.define ("caltable", caltbl+".A.cal");
1135 0 : simpar.define ("mode", mode);
1136 0 : simpar.define ("combine", ""); // SPW,FIELD, etc
1137 :
1138 0 : 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 0 : } else if (mode=="tsys-atm" or mode=="tsys-manual") {
1145 : // do be scaled in a minute by a Tsys-derived M below
1146 0 : simpar.define ("amplitude", Float(1.0) );
1147 : // os << "adding noise with unity amplitude" << LogIO::POST;
1148 0 : simpar.define ("mode", mode);
1149 :
1150 : } else {
1151 0 : throw(AipsError("unsupported mode "+mode+" in setnoise()"));
1152 : }
1153 :
1154 0 : Bool saveOnthefly=false;
1155 0 : if (simpar.isDefined("onthefly")) {
1156 0 : saveOnthefly=simpar.asBool("onthefly");
1157 : }
1158 0 : simpar.define("onthefly",OTF);
1159 :
1160 : // create the ANoise
1161 0 : if (!create_corrupt(simpar))
1162 0 : throw(AipsError("could not create ANoise in Simulator::setnoise"));
1163 :
1164 0 : if (mode=="tsys-atm" or mode=="tsys-manual") {
1165 :
1166 0 : simpar.define ("antefficiency" ,antefficiency );
1167 0 : simpar.define ("correfficiency" ,correfficiency );
1168 0 : simpar.define ("spillefficiency",spillefficiency);
1169 0 : simpar.define ("trx" ,trx );
1170 0 : simpar.define ("tground" ,tground );
1171 0 : simpar.define ("tcmb" ,tcmb );
1172 :
1173 0 : if (rxtype>=0) {
1174 0 : simpar.define ("rxType", rxtype);
1175 0 : 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 0 : if ( senscoeff > 0.0 ) {
1184 0 : simpar.define ("senscoeff", Float(senscoeff) );
1185 0 : os << "adding noise with the sensitivity constant of " << senscoeff << LogIO::POST;
1186 : } else {
1187 0 : simpar.define("senscoeff", Float(1./ sqrt(2.)));
1188 0 : os << "adding noise with the sensitivity constant of 1/sqrt(2)" << LogIO::POST;
1189 : }
1190 :
1191 0 : if (pwv.getValue("mm")>0.) {
1192 0 : 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 0 : simpar.define ("mean_pwv", pwv.getValue("mm"));
1196 : } else {
1197 0 : simpar.define ("mean_pwv", Double(1.));
1198 : // we want to set it, but it doesn't get used unless ATM is being used
1199 0 : if (mode=="tsys-atm") os<<"User has not set PWV, using 1mm"<<LogIO::POST;
1200 : }
1201 :
1202 0 : if (mode=="tsys-manual") {
1203 : // user can override the ATM calculated optical depth
1204 : // with tau0 to be used over the entire SPW,
1205 0 : simpar.define ("tau0" ,tau );
1206 0 : if (tatmos>10)
1207 0 : 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 0 : 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 0 : if (pground.getValue("mbar")>100.)
1223 0 : simpar.define ("pground", pground.getValue("mbar"));
1224 : else {
1225 0 : simpar.define ("pground", Double(560.));
1226 0 : os<<"User has not set ground pressure, using 560mb"<<LogIO::POST;
1227 : }
1228 0 : if (relhum>0)
1229 0 : 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 0 : if (altitude.getValue("m")>0.)
1235 0 : simpar.define ("altitude", altitude.getValue("m"));
1236 : else {
1237 0 : simpar.define ("altitude", Double(5000.));
1238 0 : os<<"User has not set site altitude, using 5000m"<<LogIO::POST;
1239 : }
1240 0 : if (waterheight.getValue("m")>100.)
1241 0 : simpar.define ("waterheight", waterheight.getValue("km"));
1242 : else {
1243 0 : simpar.define ("waterheight", Double(2.0));
1244 0 : 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 0 : simpar.define ("type", "TF NOISE");
1249 : // simpar.define ("type", "MF");
1250 : }
1251 :
1252 0 : if (strlen>1)
1253 0 : simpar.define ("caltable", caltbl+".T.cal");
1254 : //simpar.define ("caltable", caltbl+".M.cal");
1255 :
1256 0 : simpar.define("onthefly",saveOnthefly);
1257 :
1258 : // create the T
1259 0 : if (!create_corrupt(simpar))
1260 0 : throw(AipsError("could not create T in Simulator::setnoise"));
1261 : //throw(AipsError("could not create M in Simulator::setnoise"));
1262 : }
1263 :
1264 0 : } catch (AipsError x) {
1265 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1266 0 : return false;
1267 0 : }
1268 0 : return true;
1269 0 : }
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 0 : Bool Simulator::create_corrupt(Record& simpar)
1712 : {
1713 0 : LogIO os(LogOrigin("Simulator", "create_corrupt()", WHERE));
1714 0 : SolvableVisCal *svc(NULL);
1715 :
1716 : // RI todo sim::create_corrupt assert that ms has certain structure
1717 :
1718 : try {
1719 0 : makeVisSet();
1720 :
1721 0 : String upType=simpar.asString("type");
1722 0 : upType.upcase();
1723 0 : os << LogIO::NORMAL << "Creating "<< upType <<" Calibration structure for data corruption." << LogIO::POST;
1724 :
1725 0 : svc = createSolvableVisCal(upType,*vs_p);
1726 :
1727 0 : svc->setPrtlev(prtlev());
1728 :
1729 0 : Vector<Double> solTimes;
1730 0 : svc->setSimulate(*vs_p,simpar,solTimes);
1731 :
1732 : // add to the pointer block of VCs:
1733 0 : uInt napp=vc_p.nelements();
1734 0 : vc_p.resize(napp+1,false,true);
1735 0 : vc_p[napp] = (VisCal*) svc;
1736 : // svc=NULL;
1737 0 : ve_p.setapply(vc_p);
1738 :
1739 0 : } catch (AipsError x) {
1740 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
1741 0 : if (svc) delete svc;
1742 0 : throw(AipsError("Error in Simulator::create_corrupt"));
1743 : return false;
1744 0 : }
1745 0 : return true;
1746 0 : }
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 0 : Bool Simulator::corrupt() {
1864 :
1865 : // VIS-plane (only) corruption
1866 :
1867 0 : LogIO os(LogOrigin("Simulator", "corrupt()", WHERE));
1868 :
1869 : try {
1870 :
1871 0 : ms_p->lock();
1872 0 : 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 0 : Block<Int> columns;
1881 : // include scan iteration
1882 0 : columns.resize(5);
1883 0 : columns[0]=MS::ARRAY_ID;
1884 0 : columns[1]=MS::SCAN_NUMBER;
1885 0 : columns[2]=MS::FIELD_ID;
1886 0 : columns[3]=MS::DATA_DESC_ID;
1887 0 : columns[4]=MS::TIME;
1888 :
1889 0 : if(vs_p) {
1890 0 : delete vs_p; vs_p=0;
1891 : }
1892 0 : Matrix<Int> noselection;
1893 0 : if(mssel_p) {
1894 0 : vs_p = new VisSet(*mssel_p,columns,noselection);
1895 : }
1896 : else {
1897 0 : vs_p = new VisSet(*ms_p,columns,noselection);
1898 : }
1899 0 : 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 0 : VisIter& vi=vs_p->iter();
1906 0 : 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 0 : 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 0 : ve_p.setPivot(VisCal::B);
1915 :
1916 : // Apply
1917 0 : if (vc_p.nelements()>0) {
1918 : os << LogIO::NORMAL << "Doing visibility corruption."
1919 0 : << LogIO::POST;
1920 0 : 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 0 : os << vc_p[i]->siminfo() << "spwok = " << vc_p[i]->spwOK();
1924 0 : if (vc_p[i]->type() >= ve_p.pivot())
1925 0 : os << " in corrupt mode." << endl << LogIO::POST;
1926 : else
1927 0 : os << " in correct mode." << endl << LogIO::POST;
1928 : }
1929 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1930 : // Only if
1931 0 : if (ve_p.spwOK(vi.spectralWindow())) {
1932 0 : for (vi.origin(); vi.more(); vi++) {
1933 :
1934 : // corrupt model to pivot, correct data up to pivot
1935 0 : ve_p.collapseForSim(vb);
1936 :
1937 : // Deposit corrupted visibilities into DATA
1938 : // vi.setVis(vb.modelVisCube(), VisibilityIterator::Observed);
1939 0 : 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 0 : vi.setVis(vb.visCube(), VisibilityIterator::Corrected);
1946 :
1947 : // RI TODO is this 100% right?
1948 0 : 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 0 : 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 0 : 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 0 : ms_p->relinquishAutoLocks();
1983 0 : ms_p->unlock();
1984 0 : if(mssel_p) mssel_p->unlock();
1985 :
1986 0 : } catch (std::exception& x) {
1987 0 : ms_p->relinquishAutoLocks();
1988 0 : ms_p->unlock();
1989 0 : if(mssel_p) mssel_p->unlock();
1990 0 : throw(AipsError(x.what()));
1991 : return false;
1992 0 : }
1993 0 : return true;
1994 0 : }
1995 :
1996 :
1997 :
1998 :
1999 :
2000 :
2001 :
2002 :
2003 :
2004 :
2005 :
2006 :
2007 :
2008 :
2009 :
2010 :
2011 :
2012 :
2013 :
2014 :
2015 2 : Bool Simulator::observe(const String& sourcename,
2016 : const String& spwname,
2017 : const Quantity& startTime,
2018 : const Quantity& stopTime,
2019 : const Bool add_observation,
2020 : const Bool state_sig,
2021 : const Bool state_ref,
2022 : const double& state_cal,
2023 : const double& state_load,
2024 : const unsigned int state_sub_scan,
2025 : const String& state_obs_mode,
2026 : const String& observername,
2027 : const String& projectname)
2028 : {
2029 4 : LogIO os(LogOrigin("Simulator", "observe()", WHERE));
2030 :
2031 :
2032 : try {
2033 :
2034 2 : if(!feedsHaveBeenSet && !feedsInitialized) {
2035 0 : os << "Feeds have not been set - using default " << feedMode_p << LogIO::WARN;
2036 0 : sim_p->initFeeds(feedMode_p);
2037 0 : feedsInitialized = true;
2038 : }
2039 2 : if(!timesHaveBeenSet_p) {
2040 : os << "Times have not been set - using defaults " << endl
2041 : << " Times will be interpreted as hour angles for first source"
2042 0 : << LogIO::WARN;
2043 : }
2044 :
2045 2 : sim_p->observe(sourcename, spwname, startTime, stopTime,
2046 : add_observation, state_sig, state_ref, state_cal,state_load,state_sub_scan,state_obs_mode,observername,projectname);
2047 :
2048 :
2049 2 : if(ms_p) delete ms_p; ms_p=0;
2050 2 : if(mssel_p) delete mssel_p; mssel_p=0;
2051 4 : ms_p = new MeasurementSet(msname_p,
2052 2 : TableLock(TableLock::AutoNoReadLocking),
2053 2 : Table::Update);
2054 :
2055 2 : ms_p->flush();
2056 2 : ms_p->unlock();
2057 :
2058 0 : } catch (AipsError x) {
2059 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2060 0 : return false;
2061 0 : }
2062 2 : return true;
2063 2 : }
2064 :
2065 :
2066 :
2067 :
2068 0 : 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 0 : LogIO os(LogOrigin("Simulator", "observemany()", WHERE));
2084 :
2085 :
2086 : try {
2087 :
2088 0 : 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 0 : 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 0 : 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 0 : if(ms_p) delete ms_p; ms_p=0;
2103 0 : if(mssel_p) delete mssel_p; mssel_p=0;
2104 0 : ms_p = new MeasurementSet(msname_p,
2105 0 : TableLock(TableLock::AutoNoReadLocking),
2106 0 : Table::Update);
2107 :
2108 0 : ms_p->flush();
2109 0 : 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 0 : return true;
2116 0 : }
2117 :
2118 :
2119 :
2120 :
2121 :
2122 :
2123 0 : Bool Simulator::predict(const Vector<String>& modelImage,
2124 : const String& compList,
2125 : const Bool incremental) {
2126 :
2127 0 : 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 0 : " and componentList: " << compList << LogIO::POST;
2137 0 : if (incremental) {
2138 0 : os << "The data column will be incremented" << LogIO::POST;
2139 : } else {
2140 0 : os << "The data column will be replaced" << LogIO::POST;
2141 : }
2142 0 : if(!ms_p) {
2143 : os << "MeasurementSet pointer is null : logic problem!"
2144 0 : << LogIO::EXCEPTION;
2145 : }
2146 :
2147 0 : bool heterogeneous=False;
2148 0 : for (uInt i=1;i<diam_p.nelements();i++)
2149 0 : if (diam_p(i)!=diam_p(0)) heterogeneous=True;
2150 0 : 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 0 : ms_p->lock();
2160 0 : if(mssel_p) mssel_p->lock();
2161 0 : if (!createSkyEquation( modelImage, compList)) {
2162 0 : os << LogIO::SEVERE << "Failed to create SkyEquation" << LogIO::POST;
2163 0 : return false;
2164 : }
2165 0 : if (incremental) {
2166 0 : se_p->predict(false,MS::MODEL_DATA);
2167 : } else {
2168 0 : se_p->predict(false,MS::DATA); //20091030 RI changed SE::predict to use DATA
2169 : }
2170 0 : destroySkyEquation();
2171 :
2172 : // Copy the predicted visibilities over to the observed and
2173 : // the corrected data columns
2174 0 : makeVisSet();
2175 :
2176 0 : VisIter& vi = vs_p->iter();
2177 0 : VisBuffer vb(vi);
2178 0 : vi.origin();
2179 0 : vi.originChunks();
2180 :
2181 : //os << "Copying predicted visibilities from MODEL_DATA to DATA and CORRECTED_DATA" << LogIO::POST;
2182 :
2183 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()){
2184 0 : for (vi.origin(); vi.more(); vi++) {
2185 : // vb.setVisCube(vb.modelVisCube());
2186 :
2187 0 : 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 0 : vi.setVis(vb.visCube(),VisibilityIterator::Corrected);
2200 : }
2201 0 : vb.setModelVisCube(Complex(1.0,0.0));
2202 0 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
2203 : }
2204 : }
2205 0 : ms_p->unlock();
2206 0 : if(mssel_p) mssel_p->lock();
2207 :
2208 0 : } catch (AipsError x) {
2209 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2210 0 : ms_p->unlock();
2211 0 : if(mssel_p) mssel_p->lock();
2212 0 : return false;
2213 0 : }
2214 0 : return true;
2215 0 : }
2216 :
2217 :
2218 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2219 : // copied from SynthesisImager
2220 0 : void Simulator::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2221 : {
2222 0 : LogIO os(LogOrigin("Simulator", "getVPRecord",WHERE));
2223 :
2224 0 : VPManager *vpman=VPManager::Instance();
2225 0 : if ((doDefaultVP_p)||(vpTableStr_p != String(""))) {
2226 0 : if (doDefaultVP_p) {
2227 0 : os << "Using default Voltage Patterns from the VPManager" << LogIO::POST;
2228 0 : vpman->reset();
2229 : } else {
2230 0 : os << "Loading Voltage Pattern information from " << vpTableStr_p << LogIO::POST;
2231 0 : vpman->loadfromtable(vpTableStr_p );
2232 : }
2233 0 : 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 0 : os << "Using Voltage Patterns from the VPManager" << LogIO::POST;
2239 : }
2240 :
2241 : // PBMath::CommonPB kpb;
2242 0 : PBMath::enumerateCommonPB(telescop, kpb);
2243 : // Record rec;
2244 0 : vpman->getvp(rec, telescop);
2245 :
2246 : //ostringstream oss; rec.print(oss);
2247 : //os << "Using VP model : " << oss.str() << LogIO::POST;
2248 :
2249 0 : 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 0 : }// get VPRecord
2263 :
2264 :
2265 :
2266 0 : Bool Simulator::createSkyEquation(const Vector<String>& image,
2267 : const String complist)
2268 : {
2269 :
2270 0 : LogIO os(LogOrigin("Simulator", "createSkyEquation()", WHERE));
2271 :
2272 : try {
2273 0 : if(sm_p==0) {
2274 0 : sm_p = new CleanImageSkyModel();
2275 : }
2276 0 : AlwaysAssert(sm_p, AipsError);
2277 :
2278 : // Add the componentlist
2279 0 : if(complist!="") {
2280 0 : if(!Table::isReadable(complist)) {
2281 : os << LogIO::SEVERE << "ComponentList " << complist
2282 0 : << " not readable" << LogIO::POST;
2283 0 : return false;
2284 : }
2285 0 : componentList_p=new ComponentList(complist, true);
2286 0 : if(componentList_p==0) {
2287 : os << LogIO::SEVERE << "Cannot create ComponentList from " << complist
2288 0 : << LogIO::POST;
2289 0 : return false;
2290 : }
2291 0 : 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 0 : componentList_p=0;
2298 : }
2299 :
2300 0 : nmodels_p = image.nelements();
2301 0 : if (nmodels_p == 1 && image(0) == "") nmodels_p = 0;
2302 :
2303 0 : int models_found=0;
2304 0 : if (nmodels_p > 0) {
2305 0 : images_p.resize(nmodels_p);
2306 :
2307 0 : for (Int model=0;model<Int(nmodels_p);model++) {
2308 0 : 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 0 : if(!Table::isReadable(image(model))) {
2314 0 : os << LogIO::SEVERE << image(model) << " is unreadable" << LogIO::POST;
2315 : } else {
2316 0 : images_p[model]=0;
2317 : os << "Opening model " << model << " named "
2318 0 : << image(model) << LogIO::POST;
2319 0 : images_p[model]=new PagedImage<Float>(image(model));
2320 0 : AlwaysAssert(images_p[model], AipsError);
2321 : // RI TODO is this a logic problem with more than one source??
2322 : // Add distance
2323 0 : 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 0 : CoordinateSystem cs = images_p[model]->coordinates();
2333 0 : String errorMsg;
2334 0 : cs.setSpectralConversion(errorMsg,MFrequency::showType(MFrequency::LSRK));
2335 0 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
2336 0 : AlwaysAssert(spectralIndex>=0, AipsError);
2337 0 : SpectralCoordinate spectralCoord=cs.spectralCoordinate(spectralIndex);
2338 0 : Vector<String> units(1); units = "Hz";
2339 : // doesn't work: cs.spectralCoordinate(spectralIndex).setWorldAxisUnits(units);
2340 0 : spectralCoord.setWorldAxisUnits(units);
2341 : // put spectralCoord back into cs
2342 0 : cs.replaceCoordinate(spectralCoord,spectralIndex);
2343 : // and cs into image
2344 0 : images_p[model]->setCoordinateInfo(cs);
2345 :
2346 :
2347 0 : 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 0 : models_found++;
2352 0 : }
2353 : }
2354 : }
2355 : }
2356 :
2357 0 : 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 0 : if(vs_p) {
2365 0 : delete vs_p; vs_p=0;
2366 : }
2367 0 : makeVisSet();
2368 :
2369 0 : cft_p = new SimpleComponentFTMachine();
2370 :
2371 0 : MeasurementSet *ams=0;
2372 :
2373 0 : if(mssel_p) {
2374 0 : ams=mssel_p;
2375 : }
2376 : else {
2377 0 : ams=ms_p;
2378 : }
2379 : // 2016 from SynthesisImager:
2380 0 : MSColumns msc(*ams);
2381 0 : String telescop=msc.observation().telescopeName()(0);
2382 : PBMath::CommonPB kpb;
2383 0 : Record rec;
2384 0 : getVPRecord( rec, kpb, telescop );
2385 :
2386 0 : if((ftmachine_p=="sd")||(ftmachine_p=="both")||(ftmachine_p=="mosaic")) {
2387 0 : if(!gvp_p) {
2388 0 : if (rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2389 0 : String commonPBName;
2390 0 : PBMath::nameCommonPB(kpb, commonPBName);
2391 0 : os << "Using common PB "<<commonPBName<<" for beam calculation for telescope " << telescop << LogIO::POST;
2392 0 : gvp_p=new VPSkyJones(msc, true, parAngleInc_p, squintType_p, skyPosThreshold_p);
2393 0 : }
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 0 : if(ftmachine_p=="sd") {
2404 0 : os << "Performing Single dish gridding " << LogIO::POST;
2405 0 : if(gridfunction_p=="pb") {
2406 0 : 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 0 : 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 0 : vi.setRowBlocking(100);
2454 : }
2455 :
2456 : else {
2457 0 : os << "Synthesis gridding " << LogIO::POST;
2458 : // Now make the FTMachine
2459 : // if(wprojPlanes_p>1) {
2460 0 : 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 0 : 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 0 : os << "Fourier transforms will use image centers as tangent points" << LogIO::POST;
2535 0 : ft_p = new GridFT(cache_p/2, tile_p, gridfunction_p, mLocation_p, padding_p);
2536 : }
2537 : }
2538 0 : 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 0 : 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 0 : if(doVP_p) {
2550 0 : 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 0 : 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 0 : vp_p->summary();
2564 0 : se_p->setSkyJones(*vp_p);
2565 : } else {
2566 0 : vp_p=0;
2567 : }
2568 0 : } catch (AipsError x) {
2569 0 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2570 0 : ms_p->unlock();
2571 0 : if(mssel_p) mssel_p->lock();
2572 0 : return false;
2573 0 : }
2574 0 : return true;
2575 0 : };
2576 :
2577 0 : void Simulator::destroySkyEquation()
2578 : {
2579 0 : if(se_p) delete se_p; se_p=0;
2580 0 : if(sm_p) delete sm_p; sm_p=0;
2581 0 : if(vp_p) delete vp_p; vp_p=0;
2582 0 : if(componentList_p) delete componentList_p; componentList_p=0;
2583 :
2584 0 : for (Int model=0;model<Int(nmodels_p);model++) {
2585 0 : if(images_p[model]) delete images_p[model]; images_p[model]=0;
2586 : }
2587 0 : };
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 2 : Bool Simulator::setauto(const Double autocorrwt)
2620 : {
2621 :
2622 4 : LogIO os(LogOrigin("Simulator", "setauto()", WHERE));
2623 :
2624 : try {
2625 :
2626 2 : sim_p->setAutoCorrelationWt(autocorrwt);
2627 :
2628 : } catch (AipsError x) {
2629 : os << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2630 : return false;
2631 : }
2632 2 : return true;
2633 2 : }
2634 :
2635 0 : void Simulator::makeVisSet() {
2636 :
2637 0 : if(vs_p) return;
2638 :
2639 0 : Block<int> sort(0);
2640 0 : sort.resize(5);
2641 0 : sort[0] = MS::FIELD_ID;
2642 0 : sort[1] = MS::FEED1;
2643 0 : sort[2] = MS::ARRAY_ID;
2644 0 : sort[3] = MS::DATA_DESC_ID;
2645 0 : sort[4] = MS::TIME;
2646 0 : Matrix<Int> noselection;
2647 0 : if(mssel_p) {
2648 0 : 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 0 : AlwaysAssert(vs_p, AipsError);
2657 :
2658 0 : }
2659 :
2660 :
2661 :
2662 0 : 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 0 : LogIO os(LogOrigin("Simulator", "setoptions()", WHERE));
2668 :
2669 0 : os << "Setting processing options" << LogIO::POST;
2670 :
2671 0 : sim_p->setMaxData(maxData*1024.0*1024.0);
2672 :
2673 0 : ftmachine_p=downcase(ftmachine);
2674 0 : if(cache>0) cache_p=cache;
2675 0 : if(tile>0) tile_p=tile;
2676 0 : gridfunction_p=downcase(gridfunction);
2677 0 : mLocation_p=mLocation;
2678 0 : if(padding>=1.0) {
2679 0 : padding_p=padding;
2680 : }
2681 0 : facets_p=facets;
2682 0 : wprojPlanes_p = wprojPlanes;
2683 : // Destroy the FTMachine
2684 0 : if(ft_p) {delete ft_p; ft_p=0;}
2685 0 : if(cft_p) {delete cft_p; cft_p=0;}
2686 :
2687 0 : return true;
2688 0 : }
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 0 : Bool Simulator::setdata(const Vector<Int>& spectralwindowids,
2717 : const Vector<Int>& fieldids,
2718 : const String& msSelect)
2719 :
2720 : {
2721 :
2722 :
2723 0 : LogIO os(LogOrigin("Simulator", "setdata()", WHERE));
2724 :
2725 0 : 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 0 : os << "Selecting data" << LogIO::POST;
2734 :
2735 : // Map the selected spectral window ids to data description ids
2736 0 : MSDataDescColumns dataDescCol(ms_p->dataDescription());
2737 0 : Vector<Int> ddSpwIds=dataDescCol.spectralWindowId().getColumn();
2738 :
2739 0 : Vector<Int> datadescids(0);
2740 0 : for (uInt row=0; row<ddSpwIds.nelements(); row++) {
2741 0 : Bool found=false;
2742 0 : for (uInt j=0; j<spectralwindowids.nelements(); j++) {
2743 0 : if (ddSpwIds(row)==spectralwindowids(j)) found=true;
2744 : };
2745 0 : if (found) {
2746 0 : datadescids.resize(datadescids.nelements()+1,true);
2747 0 : datadescids(datadescids.nelements()-1)=row;
2748 : };
2749 : };
2750 :
2751 0 : if(vs_p) delete vs_p; vs_p=0;
2752 0 : 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 0 : if(fieldids.nelements()>0||datadescids.nelements()>0) {
2758 0 : os << "Performing selection on MeasurementSet" << LogIO::POST;
2759 0 : Table& original=*ms_p;
2760 :
2761 : // Now we make a condition to do the old FIELD_ID, SPECTRAL_WINDOW_ID
2762 : // selection
2763 0 : TableExprNode condition;
2764 0 : String colf=MS::columnName(MS::FIELD_ID);
2765 0 : String cols=MS::columnName(MS::DATA_DESC_ID);
2766 0 : if(fieldids.nelements()>0&&datadescids.nelements()>0){
2767 0 : condition=original.col(colf).in(fieldids)&&original.col(cols).in(datadescids);
2768 0 : os << "Selecting on field and spectral window ids" << LogIO::POST;
2769 : }
2770 0 : else if(datadescids.nelements()>0) {
2771 0 : condition=original.col(cols).in(datadescids);
2772 0 : 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 0 : mssel_p = new MeasurementSet(original(condition));
2781 :
2782 : //AlwaysAssert(mssel_p, AipsError);
2783 : //mssel_p->rename(msname_p+"/SELECTED_TABLE", Table::Scratch);
2784 0 : 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 0 : mssel_p->flush();
2793 : }
2794 :
2795 0 : }
2796 : else {
2797 0 : mssel_p=new MeasurementSet(*ms_p);
2798 : }
2799 : {
2800 0 : Int len = msSelect.length();
2801 0 : Int nspace = msSelect.freq (' ');
2802 0 : Bool nullSelect=(msSelect.empty() || nspace==len);
2803 0 : 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 0 : os << "No selection string given" << LogIO::POST;
2827 : }
2828 :
2829 0 : if(mssel_p->nrow()!=ms_p->nrow()) {
2830 0 : os << "By selection " << ms_p->nrow() << " rows are reduced to "
2831 0 : << mssel_p->nrow() << LogIO::POST;
2832 : }
2833 : else {
2834 0 : os << "Selection did not drop any rows" << LogIO::POST;
2835 : }
2836 : }
2837 :
2838 : // Now create the VisSet
2839 0 : if(vs_p) delete vs_p; vs_p=0;
2840 0 : makeVisSet();
2841 : //Now assign the source directions to something selected or sensible
2842 : {
2843 : //Int fieldsel=0;
2844 0 : if(fieldids.nelements() >0) {
2845 : //fieldsel=fieldids(0);
2846 : // RI TODO does sim:setdata need this?
2847 0 : nField=fieldids.nelements();
2848 0 : 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 0 : (vs_p->iter()).msColumns().field().name().get(i,sourceName_p[i]);
2852 0 : sourceDirection_p[i]=(vs_p->iter()).msColumns().field().phaseDirMeas(i);
2853 0 : (vs_p->iter()).msColumns().field().code().get(i,calCode_p[i]);
2854 : }
2855 : }
2856 : }
2857 0 : return true;
2858 0 : } 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 0 : }
2865 :
2866 : } // end namespace casa
|
