Line data Source code
1 : //# SimpleSimVi2.cc: Rudimentary data simulator--implementation
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 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: VisibilityIterator2.h,v 19.14 2006/02/28 04:48:58 mvoronko Exp $
27 :
28 : #include <msvis/MSVis/SimpleSimVi2.h>
29 : #include <casacore/measures/Measures/MFrequency.h>
30 : #include <casacore/measures/Measures/MEpoch.h>
31 : #include <casacore/casa/Arrays.h>
32 : #include <casacore/casa/BasicMath/Random.h>
33 : #include <casacore/casa/Quanta/MVTime.h>
34 : #include <casacore/casa/Containers/Record.h>
35 : #include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
36 : #include <casacore/tables/Tables/SetupNewTab.h>
37 :
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : namespace vi {
42 :
43 22 : SimpleSimVi2Parameters::SimpleSimVi2Parameters() :
44 22 : nField_(1),
45 22 : nScan_(1),
46 22 : nSpw_(1),
47 22 : nAnt_(4),
48 22 : nCorr_(4),
49 22 : nTimePerField_(Vector<Int>(1,1)),
50 22 : nChan_(Vector<Int>(1,1)),
51 22 : date0_("2016/01/06/00:00:00."),
52 22 : dt_(1.0),
53 22 : refFreq_(Vector<Double>(1,100.0e9)),
54 22 : df_(Vector<Double>(1,1.0e6)), // 1 MHz chans
55 22 : doNoise_(false),
56 22 : stokes_(Matrix<Float>(1,1,1.0)),
57 22 : gain_(Matrix<Float>(1,1,1.0)),
58 22 : tsys_(Matrix<Float>(1,1,1.0)),
59 22 : doNorm_(false),
60 22 : polBasis_("circ"),
61 22 : doAC_(false),
62 22 : c0_(Complex(0.0)),
63 22 : autoPol_(false),
64 22 : doParang_(false),
65 22 : spwScope_(ChunkScope),
66 22 : antennaScope_(RowScope)
67 : {
68 :
69 22 : Vector<Int> nTimePerField(nTimePerField_);
70 22 : Vector<Int> nChan(nChan_);
71 22 : Vector<Double> refFreq(refFreq_);
72 22 : Vector<Double> df(df_);
73 22 : Matrix<Float> stokes(stokes_);
74 22 : Matrix<Float> gain(gain_);
75 22 : Matrix<Float> tsys(tsys_);
76 :
77 :
78 : // Generic initialization
79 22 : this->initialize(nTimePerField,nChan,refFreq,df,stokes,gain,tsys);
80 :
81 22 : }
82 :
83 0 : SimpleSimVi2Parameters::SimpleSimVi2Parameters(Int nField,Int nScan,Int nSpw,Int nAnt,Int nCorr,
84 : const Vector<Int>& nTimePerField,const Vector<Int>& nChan,
85 : Complex c0,
86 : String polBasis,
87 : Bool autoPol,Bool doParang,
88 0 : Bool doAC) :
89 0 : nField_(nField),
90 0 : nScan_(nScan),
91 0 : nSpw_(nSpw),
92 0 : nAnt_(nAnt),
93 0 : nCorr_(nCorr),
94 0 : nTimePerField_(),
95 0 : nChan_(),
96 0 : date0_("2016/01/06/00:00:00."), // the rest are defaulted
97 0 : dt_(1.0),
98 0 : refFreq_(Vector<Double>(1,100.0e9)),
99 0 : df_(Vector<Double>(nSpw,1.0e6)), // 1 MHz chans
100 0 : doNoise_(False),
101 0 : stokes_(Matrix<Float>(1,1,1.0)),
102 0 : gain_(Matrix<Float>(1,1,1.0)),
103 0 : tsys_(Matrix<Float>(1,1,1.0)),
104 0 : doNorm_(False),
105 0 : polBasis_(polBasis),
106 0 : doAC_(doAC),
107 0 : c0_(c0),
108 0 : autoPol_(autoPol),
109 0 : doParang_(doParang),
110 0 : spwScope_(ChunkScope),
111 0 : antennaScope_(RowScope)
112 : {
113 :
114 0 : Vector<Double> refFreq(nSpw_,100.0e9);
115 : // Make multiple spws adjacent (not identical) in frequency
116 0 : Int uNChan=nChan.nelements();
117 0 : for (Int i=1;i<nSpw_;i++)
118 0 : refFreq[i]=refFreq[i-1]+Double(nChan[(i-1)%uNChan]*df_[i-1]);
119 :
120 : // cout << "SSV::refFreq=" << refFreq << endl;
121 :
122 0 : Vector<Double> df(df_);
123 0 : Matrix<Float> stokes(stokes_);
124 0 : Matrix<Float> gain(gain_);
125 0 : Matrix<Float> tsys(tsys_);
126 :
127 : // cout << "SSVP(simple): stokes = " << stokes << " " << " doParang_=" << boolalpha << doParang_ << endl;
128 :
129 : // Generic initialization
130 0 : this->initialize(nTimePerField,nChan,refFreq,df,stokes,gain,tsys);
131 :
132 0 : }
133 :
134 :
135 :
136 0 : SimpleSimVi2Parameters::SimpleSimVi2Parameters(Int nField,Int nScan,Int nSpw,Int nAnt,Int nCorr,
137 : const Vector<Int>& nTimePerField,const Vector<Int>& nChan,
138 : String date0, Double dt,
139 : const Vector<Double>& refFreq, const Vector<Double>& df,
140 : const Matrix<Float>& stokes,
141 : Bool doNoise,
142 : const Matrix<Float>& gain, const Matrix<Float>& tsys,
143 : Bool doNorm,
144 : String polBasis, Bool doAC,
145 : Complex c0, Bool doParang,
146 : MetadataScope spwScope,
147 0 : MetadataScope antennaScope) :
148 0 : nField_(nField),
149 0 : nScan_(nScan),
150 0 : nSpw_(nSpw),
151 0 : nAnt_(nAnt),
152 0 : nCorr_(nCorr),
153 0 : nTimePerField_(),
154 0 : nChan_(),
155 0 : date0_(date0),
156 0 : dt_(dt),
157 0 : refFreq_(),
158 0 : df_(),
159 0 : doNoise_(doNoise),
160 0 : stokes_(),
161 0 : gain_(),
162 0 : tsys_(),
163 0 : doNorm_(doNorm),
164 0 : polBasis_(polBasis),
165 0 : doAC_(doAC),
166 0 : c0_(c0),
167 0 : autoPol_(false),
168 0 : doParang_(doParang),
169 0 : spwScope_(spwScope),
170 0 : antennaScope_(antennaScope)
171 : {
172 :
173 : // Generic initialization
174 0 : this->initialize(nTimePerField,nChan,refFreq,df,stokes,gain,tsys);
175 :
176 0 : }
177 :
178 0 : SimpleSimVi2Parameters::SimpleSimVi2Parameters(const SimpleSimVi2Parameters& other) {
179 0 : *this=other;
180 0 : }
181 :
182 0 : SimpleSimVi2Parameters& SimpleSimVi2Parameters::operator=(const SimpleSimVi2Parameters& other) {
183 :
184 0 : if (this != &other) {
185 0 : nField_=other.nField_;
186 0 : nScan_=other.nScan_;
187 0 : nSpw_=other.nSpw_;
188 0 : nAnt_=other.nAnt_;
189 0 : nCorr_=other.nCorr_;
190 0 : nTimePerField_.assign(other.nTimePerField_); // NB: Array::assign() forces reshape
191 0 : nChan_.assign(other.nChan_);
192 0 : date0_=other.date0_;
193 0 : dt_=other.dt_;
194 0 : refFreq_.assign(other.refFreq_);
195 0 : df_.assign(other.df_);
196 0 : doNoise_=other.doNoise_;
197 0 : stokes_.assign(other.stokes_);
198 0 : gain_.assign(other.gain_);
199 0 : tsys_.assign(other.tsys_);
200 0 : doNorm_=other.doNorm_;
201 0 : polBasis_=other.polBasis_;
202 0 : doAC_=other.doAC_;
203 0 : c0_=other.c0_;
204 0 : autoPol_=other.autoPol_;
205 0 : doParang_=other.doParang_;
206 0 : spwScope_=other.spwScope_;
207 0 : antennaScope_=other.antennaScope_;
208 : }
209 0 : return *this;
210 :
211 : }
212 :
213 :
214 22 : SimpleSimVi2Parameters::~SimpleSimVi2Parameters() {}
215 :
216 :
217 :
218 0 : void SimpleSimVi2Parameters::summary() const {
219 :
220 0 : cout << endl << "***SimpleSimVi2Parameters Summary******************" << endl;
221 0 : cout << boolalpha;
222 0 : cout << "* nField = " << nField_ << endl;
223 0 : cout << "* nScan = " << nScan_ << endl;
224 0 : cout << "* nSpw = " << nSpw_ << endl;
225 0 : cout << "* nAnt = " << nAnt_ << endl;
226 0 : cout << "* nCorr = " << nCorr_ << endl;
227 :
228 0 : cout << "* nTimePerField = " << nTimePerField_ << endl;
229 0 : cout << "* nChan = " << nChan_ << endl;
230 :
231 0 : cout << "* date0 = " << date0_ << endl;
232 0 : cout << "* dt = " << dt_ << endl;
233 :
234 0 : cout << "* refFreq = " << refFreq_ << endl;
235 0 : cout << "* df = " << df_ << endl;
236 :
237 0 : cout << "* stokes = " << stokes_ << endl;
238 :
239 0 : cout << "* doNoise = " << doNoise_ << endl;
240 :
241 0 : cout << "* gain = " << gain_ << endl;
242 0 : cout << "* tsys = " << tsys_ << endl;
243 :
244 0 : cout << "* doNorm = " << doNorm_ << endl;
245 :
246 0 : cout << "* polBasis = " << polBasis_ << endl;
247 0 : cout << "* doAC = " << doAC_ << endl;
248 0 : cout << "* c0 = " << c0_ << endl;
249 0 : cout << "***************************************************" << endl << endl;
250 0 : }
251 :
252 : // Return frequencies for specified spw
253 0 : Vector<Double> SimpleSimVi2Parameters::freqs(Int spw) const {
254 0 : Vector<Double> f(nChan_(spw));
255 0 : indgen(f);
256 0 : f*=df_(spw);
257 0 : f+=(df_(spw)/2. + refFreq_(spw));
258 0 : return f;
259 0 : }
260 :
261 :
262 :
263 22 : void SimpleSimVi2Parameters::initialize(const Vector<Int>& nTimePerField,const Vector<Int>& nChan,
264 : const Vector<Double>& refFreq, const Vector<Double>& df,
265 : const Matrix<Float>& stokes,
266 : const Matrix<Float>& gain, const Matrix<Float>& tsys) {
267 :
268 22 : nTimePerField_.resize(nField_); // field-dep scan length
269 22 : if (nTimePerField.nelements()==1)
270 22 : nTimePerField_.set(nTimePerField(0)); // all to specified value
271 : else
272 0 : nTimePerField_=nTimePerField; // will throw if length mismatch
273 :
274 22 : nChan_.resize(nSpw_);
275 22 : if (nChan.nelements()==1)
276 22 : nChan_.set(nChan(0)); // all to specified value
277 : else
278 0 : nChan_=nChan; // will throw if length mismatch
279 :
280 22 : refFreq_.resize(nSpw_);
281 22 : if (refFreq.nelements()==1)
282 22 : refFreq_.set(refFreq(0));
283 : else
284 0 : refFreq_=refFreq; // will throw if length mismatch
285 :
286 22 : df_.resize(nSpw_);
287 22 : if (df.nelements()==1)
288 22 : df_.set(df(0)); // all to specified value
289 : else
290 0 : df_=df; // will throw if length mismatch
291 :
292 22 : stokes_.resize(4,nField_);
293 22 : if (stokes.nelements()==1) {
294 22 : stokes_.set(0.0f); // enforce unpolarized!
295 22 : stokes_(Slice(0),Slice())=stokes(0,0); // only I specified
296 :
297 : // If requested, set Q=0.04*I, U=0.03*I
298 22 : if (autoPol_) {
299 0 : stokes_(Slice(1),Slice())=Float(stokes(0,0)*0.04f);
300 0 : stokes_(Slice(2),Slice())=Float(stokes(0,0)*0.03f);
301 : }
302 : }
303 : else
304 0 : stokes_=stokes; // insist shapes match
305 :
306 22 : gain_.resize(2,nAnt_);
307 22 : if (gain.nelements()==1)
308 22 : gain_.set(gain(0,0)); // all to specified value
309 : else
310 0 : gain_=gain; // will throw if shapes mismatch
311 :
312 :
313 22 : tsys_.resize(2,nAnt_);
314 22 : if (tsys.nelements()==1)
315 22 : tsys_.set(tsys(0,0)); // all to specified value
316 : else
317 0 : tsys_=tsys; // will throw if shapes mismatch
318 :
319 22 : if(antennaScope_ == ChunkScope)
320 0 : SSVi2NotYetImplemented();
321 22 : }
322 :
323 :
324 0 : SimpleSimVi2::SimpleSimVi2 () {}
325 :
326 0 : SimpleSimVi2::SimpleSimVi2 (const SimpleSimVi2Parameters& pars)
327 :
328 : : ViImplementation2(),
329 0 : pars_(pars),
330 0 : nChunk_(0),
331 0 : nBsln_(0),
332 0 : t0_(0.0),
333 0 : wt0_(),
334 0 : vis0_(),
335 0 : iChunk_(0),
336 0 : iSubChunk_(0),
337 0 : iRow0_(0),
338 0 : iScan_(0),
339 0 : iChunkTime0_(0),
340 0 : thisScan_(1),
341 0 : thisField_(0),
342 0 : thisSpw_(0),
343 0 : lastScan_(-1),
344 0 : lastField_(-1),
345 0 : lastSpw_(-1),
346 0 : thisTime_(0.0),
347 0 : corrdef_(4,Stokes::Undefined),
348 0 : vb_(nullptr),
349 0 : phaseCenter_(),
350 0 : feedpa_()
351 : {
352 : // Derived stuff
353 :
354 0 : if(pars_.spwScope_ == ChunkScope)
355 0 : nChunk_=pars_.nScan_*pars_.nSpw_;
356 : else
357 0 : nChunk_=pars_.nScan_;
358 :
359 0 : nBsln_=pars_.nAnt_*(pars_.nAnt_+ (pars_.doAC_ ? 1 : -1))/2;
360 :
361 : // Time tbd
362 : // cout << "Using 2016/01/06/00:00:00.0 as reference time." << endl;
363 0 : t0_=4958755200.0;
364 :
365 : // Fundamental weight value is pars_.df_*dt_
366 0 : wt0_.resize(pars_.nSpw_);
367 0 : if (pars_.doNoise_) {
368 0 : convertArray(wt0_,pars_.df_); // Float <- Double
369 0 : wt0_*=Float(pars_.dt_);
370 : }
371 : else
372 0 : wt0_.set(1.0);
373 :
374 : // Fundamental vis are just stokes combos (complex)
375 0 : const Int& nCor(pars_.nCorr_);
376 0 : vis0_.resize(nCor,pars_.nField_);
377 0 : vis0_.set(Complex(0.0));
378 0 : for (Int ifld=0;ifld<pars_.nField_;++ifld) {
379 0 : if (pars_.polBasis_=="circ") {
380 0 : vis0_(0,ifld)=Complex(pars_.stokes_(0,ifld)+pars_.stokes_(3,ifld),0.0);
381 0 : if (nCor>1) {
382 0 : vis0_(nCor-1,ifld)=Complex(pars_.stokes_(0,ifld)-pars_.stokes_(3,ifld),0.0);
383 0 : if (nCor>2) {
384 0 : vis0_(1,ifld)=Complex(pars_.stokes_(1,ifld), pars_.stokes_(2,ifld));
385 0 : vis0_(2,ifld)=Complex(pars_.stokes_(1,ifld),-1.0f*pars_.stokes_(2,ifld));
386 : }
387 : }
388 : }
389 0 : else if (pars_.polBasis_=="lin") {
390 0 : vis0_(0,ifld)=Complex(pars_.stokes_(0,ifld)+pars_.stokes_(1,ifld),0.0);
391 0 : if (nCor>1) {
392 0 : vis0_(nCor-1,ifld)=Complex(pars_.stokes_(0,ifld)-pars_.stokes_(1,ifld),0.0);
393 0 : if (nCor>2) {
394 0 : vis0_(1,ifld)=Complex(pars_.stokes_(2,ifld), pars_.stokes_(3,ifld));
395 0 : vis0_(2,ifld)=Complex(pars_.stokes_(2,ifld),-1.0f*pars_.stokes_(3,ifld));
396 : }
397 : }
398 : }
399 : }
400 :
401 0 : corrdef_.resize(nCor);
402 0 : if (pars_.polBasis_=="circ") {
403 0 : corrdef_(0)=Stokes::type("RR");
404 0 : if (nCor>1) {
405 0 : corrdef_(nCor-1)=Stokes::type("LL");
406 0 : if (nCor>2) {
407 0 : corrdef_(1)=Stokes::type("RL");
408 0 : corrdef_(2)=Stokes::type("LR");
409 : }
410 : }
411 : }
412 0 : else if (pars_.polBasis_=="lin") {
413 0 : corrdef_(0)=Stokes::type("XX");
414 0 : if (nCor>1) {
415 0 : corrdef_(3)=Stokes::type("YY");
416 0 : if (nCor>2) {
417 0 : corrdef_(1)=Stokes::type("XY");
418 0 : corrdef_(2)=Stokes::type("YX");
419 : }
420 : }
421 : }
422 :
423 0 : VisBufferOptions vbopt=VbWritable;
424 0 : vb_.reset(createAttachedVisBuffer(vbopt));
425 :
426 0 : generateSubtables();
427 0 : }
428 :
429 0 : void SimpleSimVi2::generateSubtables()
430 : {
431 : // Generating Antenna Subtable
432 0 : TableDesc antennaTD = MSAntenna::requiredTableDesc();
433 0 : SetupNewTable antennaSetup("antennaSubtable", antennaTD, Table::New);
434 0 : antennaSubTable_p = Table(antennaSetup, Table::Memory, pars_.nAnt_, true);
435 0 : antennaSubTablecols_p.reset(new MSAntennaColumns(antennaSubTable_p));
436 :
437 : // Generating SPW Subtable
438 0 : TableDesc spwTD = MSSpectralWindow::requiredTableDesc();
439 0 : SetupNewTable spwSetup("spwSubtable", spwTD, Table::New);
440 0 : spwSubTable_p = Table(spwSetup, Table::Memory, pars_.nSpw_, true);
441 0 : spwSubTablecols_p.reset(new MSSpWindowColumns(spwSubTable_p));
442 0 : auto numChanCol = spwSubTablecols_p->numChan();
443 0 : numChanCol.putColumn(pars_.nChan_);
444 0 : auto refFreqCol = spwSubTablecols_p->refFrequency();
445 0 : refFreqCol.putColumn(pars_.refFreq_);
446 0 : auto chanFreqCol = spwSubTablecols_p->chanFreq();
447 0 : for (Int i=0; i < pars_.nSpw_; i++)
448 0 : chanFreqCol.put(i, pars_.freqs(i));
449 0 : auto chanWidthCol = spwSubTablecols_p->chanWidth();
450 0 : for (Int i=0; i < pars_.nSpw_; i++)
451 0 : chanWidthCol.put(i, Vector<double>(pars_.nChan_(i), pars_.df_(i)));
452 :
453 : // Generating DD Subtable. There is only one polarizations,
454 : // therefore number of DDs = number of SPWs.
455 0 : TableDesc ddTD = MSDataDescription::requiredTableDesc();
456 0 : SetupNewTable ddSetup("ddSubtable", ddTD, Table::New);
457 0 : ddSubTable_p = Table(ddSetup, Table::Memory, pars_.nSpw_, true);
458 0 : ddSubTablecols_p.reset(new MSDataDescColumns(ddSubTable_p));
459 0 : auto spwCol = ddSubTablecols_p->spectralWindowId();
460 0 : for (int i=0; i < pars_.nSpw_; i++)
461 0 : spwCol.put(i, i);
462 :
463 : // Generating polarization Subtable. There is only one polarizations,
464 0 : TableDesc polTD = MSPolarization::requiredTableDesc();
465 0 : SetupNewTable polSetup("polSubtable", polTD, Table::New);
466 0 : polSubTable_p = Table(polSetup, Table::Memory, 1, true);
467 0 : polSubTablecols_p.reset(new MSPolarizationColumns(polSubTable_p));
468 :
469 0 : }
470 :
471 : // Destructor
472 0 : SimpleSimVi2::~SimpleSimVi2 () { /*cout << " ~SSVi2: " << this << endl;*/ }
473 :
474 :
475 : // +==================================+
476 : // | |
477 : // | Iteration Control and Monitoring |
478 : // | |
479 : // +==================================+
480 :
481 : // Methods to control chunk iterator
482 :
483 0 : void SimpleSimVi2::originChunks (Bool)
484 : {
485 : // Initialize global indices
486 0 : iChunk_=0;
487 0 : thisField_=0;
488 0 : if(pars_.spwScope_ == ChunkScope)
489 0 : thisSpw_=0;
490 :
491 : // First Scan
492 0 : thisScan_=1;
493 :
494 : // Initialize time
495 0 : iChunkTime0_=t0_+pars_.dt_/2.;
496 0 : thisTime_=iChunkTime0_;
497 :
498 0 : if(pars_.spwScope_ == ChunkScope)
499 0 : iRow0_=-nBsln_;
500 : else
501 0 : iRow0_=-pars_.nSpw_ * nBsln_;
502 :
503 0 : }
504 :
505 :
506 0 : Bool SimpleSimVi2::moreChunks () const
507 : {
508 : // if there are more chunks...
509 0 : return iChunk_<nChunk_;
510 : }
511 :
512 0 : void SimpleSimVi2::nextChunk ()
513 : {
514 : // Remember last chunk's indices
515 0 : lastScan_=thisScan_;
516 0 : lastField_=thisField_;
517 0 : if(pars_.spwScope_ == ChunkScope)
518 0 : lastSpw_=thisSpw_;
519 :
520 : // Increment chunk counter
521 0 : ++iChunk_;
522 :
523 : // If each chunk contains a unique SPW, then a new scan happens
524 : // each pars_.nSpw_ chunks
525 : // If eah chunk contains all SPWs, then a new scan happens each chunk
526 0 : if(pars_.spwScope_ == ChunkScope)
527 0 : iScan_ = iChunk_/pars_.nSpw_;
528 : else
529 0 : iScan_ = iChunk_;
530 : // 1-based
531 0 : thisScan_ = 1+ iScan_;
532 :
533 : // Each scan is new field
534 0 : thisField_ = iScan_%pars_.nField_;
535 :
536 : // Each chunk is new spw if each chunk contains a unique SPW
537 0 : if(pars_.spwScope_ == ChunkScope)
538 0 : thisSpw_ = iChunk_%pars_.nSpw_;
539 :
540 : // Increment chunk time if new scan
541 : // (spws have been exhausted on previous scan)
542 0 : if (thisScan_!=lastScan_) iChunkTime0_=thisTime_ + pars_.dt_;
543 :
544 : // Ensure subchunks initialized
545 : // this->origin();
546 :
547 0 : }
548 :
549 0 : void SimpleSimVi2::result(casacore::Record& res) const
550 : {
551 0 : res.define("comment", "SimpleSimVi2");
552 0 : }
553 :
554 : // Methods to control and monitor subchunk iteration
555 :
556 0 : void SimpleSimVi2::origin ()
557 : {
558 0 : rownr_t spwSubchunkFactor = 1, antennaSubchunkFactor = 1;
559 0 : if(pars_.spwScope_ == SubchunkScope)
560 0 : spwSubchunkFactor = pars_.nSpw_;
561 0 : if(pars_.antennaScope_ == SubchunkScope)
562 0 : antennaSubchunkFactor = nBsln_;
563 :
564 0 : nSubchunk_ = spwSubchunkFactor * antennaSubchunkFactor * pars_.nTimePerField_(thisField_);
565 :
566 : // First subchunk this chunk
567 0 : iSubChunk_=0;
568 :
569 : // time is first time of the chunk
570 0 : thisTime_=iChunkTime0_;
571 :
572 0 : rownr_t spwRowFactor = 1, antennaRowFactor = 1;
573 : // row counter
574 0 : if(pars_.spwScope_ == RowScope)
575 0 : spwRowFactor = pars_.nSpw_;
576 0 : if(pars_.antennaScope_ == RowScope)
577 0 : antennaRowFactor = nBsln_;
578 :
579 0 : iRow0_ += spwRowFactor * antennaRowFactor;
580 :
581 : // Start with SPW=0 if scope is Subchunk.
582 0 : if(pars_.spwScope_ == SubchunkScope || pars_.spwScope_ == RowScope)
583 0 : thisSpw_ = 0;
584 :
585 0 : if(pars_.antennaScope_ == SubchunkScope)
586 : {
587 0 : thisAntenna1_ = 0;
588 0 : thisAntenna2_ = (pars_.doAC_ ? 0 : 1);
589 : }
590 :
591 : // Keep VB sync'd
592 0 : this->configureNewSubchunk();
593 0 : }
594 :
595 0 : Bool SimpleSimVi2::more () const
596 : {
597 : // true if still more subchunks for this scan's field
598 0 : return (iSubChunk_<nSubchunk_);
599 : }
600 :
601 0 : void SimpleSimVi2::next () {
602 : // Advance counter and time
603 0 : ++iSubChunk_;
604 :
605 0 : if(iSubChunk_ < nSubchunk_)
606 : {
607 : // Change SPW once all times and baselines have been served
608 0 : if(pars_.spwScope_ == SubchunkScope)
609 : {
610 0 : if(pars_.antennaScope_ == SubchunkScope)
611 0 : thisSpw_ = iSubChunk_ / (nBsln_ / pars_.nTimePerField_(thisField_));
612 : else
613 0 : thisSpw_ = iSubChunk_ / pars_.nTimePerField_(thisField_);
614 : }
615 :
616 0 : if(pars_.antennaScope_ == SubchunkScope)
617 : {
618 0 : thisAntenna2_++;
619 0 : if(thisAntenna2_ == pars_.nAnt_)
620 : {
621 0 : thisAntenna1_++;
622 0 : thisAntenna2_ = (pars_.doAC_ ? thisAntenna1_ : thisAntenna1_ + 1);
623 : }
624 : // All baselines have been served, start again
625 0 : if(! (iSubChunk_ % nBsln_) )
626 : {
627 0 : thisAntenna1_ = 0;
628 0 : thisAntenna2_ = (pars_.doAC_ ? 0 : 1);
629 0 : thisTime_+=pars_.dt_;
630 : }
631 : }
632 : else
633 0 : thisTime_+=pars_.dt_;
634 :
635 : // Keep VB sync'd
636 0 : this->configureNewSubchunk();
637 : }
638 0 : }
639 :
640 0 : Subchunk SimpleSimVi2::getSubchunkId () const { return Subchunk(iChunk_,iSubChunk_);}
641 :
642 : // Return the time interval (in seconds) used for iteration.
643 : // This is not the same as the INTERVAL column. Setting the
644 : // the interval requires calling origin chunks before performing
645 : // further iterator.
646 :
647 : // Select the channels to be returned. Requires calling originChunks before
648 : // performing additional iteration.
649 :
650 0 : void SimpleSimVi2::setFrequencySelections (const FrequencySelections &)
651 : {
652 0 : SSVi2NotYetImplemented()
653 : }
654 :
655 0 : Bool SimpleSimVi2::existsColumn (VisBufferComponent2 id) const
656 : {
657 : Bool result;
658 0 : switch (id){
659 :
660 0 : case VisBufferComponent2::VisibilityCubeFloat:
661 0 : result = False;
662 0 : break;
663 :
664 0 : case VisBufferComponent2::VisibilityCorrected:
665 : case VisBufferComponent2::VisibilityCubeCorrected:
666 : case VisBufferComponent2::VisibilityModel:
667 : case VisBufferComponent2::VisibilityCubeModel:
668 : case VisBufferComponent2::VisibilityObserved:
669 : case VisBufferComponent2::VisibilityCubeObserved:
670 : case VisBufferComponent2::WeightSpectrum:
671 : case VisBufferComponent2::SigmaSpectrum:
672 0 : result = True;
673 0 : break;
674 :
675 0 : default:
676 0 : result = True; // required columns
677 0 : break;
678 : }
679 :
680 0 : return result;
681 :
682 :
683 : }
684 :
685 0 : casacore::rownr_t SimpleSimVi2::nRows () const
686 : {
687 0 : return nRows_;
688 : };
689 :
690 0 : casacore::rownr_t SimpleSimVi2::nShapes () const
691 : {
692 0 : return nShapes_;
693 : };
694 :
695 0 : const casacore::Vector<casacore::rownr_t>& SimpleSimVi2::nRowsPerShape () const
696 : {
697 0 : return nRowsPerShape_;
698 : }
699 :
700 0 : const casacore::Vector<casacore::Int>& SimpleSimVi2::nChannelsPerShape () const
701 : {
702 0 : return nChannPerShape_;
703 : }
704 :
705 0 : const casacore::Vector<casacore::Int>& SimpleSimVi2::nCorrelationsPerShape () const
706 : {
707 0 : return nCorrsPerShape_;
708 : }
709 :
710 : // Return the row ids as from the original root table. This is useful
711 : // to find correspondance between a given row in this iteration to the
712 : // original ms row
713 0 : void SimpleSimVi2::getRowIds (Vector<rownr_t> &) const {
714 0 : SSVi2NotYetImplemented()
715 : /*
716 : rowids.resize(nRows());
717 : indgen(rowids);
718 : rowids+=iRow0; // offset to current iteration
719 : */
720 : }
721 :
722 : /*
723 : VisBuffer2 * SimpleSimVi2::getVisBuffer (const VisibilityIterator2 * vi)
724 : {
725 : ThrowIf (vb_ == nullptr, "VI Implementation has not VisBuffer.");
726 : vb_->associateWithVi2 (vi);
727 : return vb_;
728 : }
729 : */
730 :
731 0 : VisBuffer2 * SimpleSimVi2::getVisBuffer () const { return vb_.get(); }
732 :
733 : // +=========================+
734 : // | |
735 : // | Subchunk Data Accessors |
736 : // | |
737 : // +=========================+
738 :
739 : // Return info
740 0 : void SimpleSimVi2::antenna1 (Vector<Int> & ant1) const {
741 0 : ant1.resize(nRows());
742 0 : if(pars_.antennaScope_ == RowScope)
743 : {
744 0 : Int k=0;
745 0 : for (Int i=0;i<(pars_.doAC_ ? pars_.nAnt_ : pars_.nAnt_-1);++i) {
746 0 : for (Int j=(pars_.doAC_ ? i : i+1);j<pars_.nAnt_;++j) {
747 0 : ant1[k]=i;
748 0 : ++k;
749 : }
750 : }
751 : }
752 : else
753 0 : ant1 = thisAntenna1_;
754 0 : }
755 0 : void SimpleSimVi2::antenna2 (Vector<Int> & ant2) const {
756 0 : ant2.resize(nRows());
757 0 : if(pars_.antennaScope_ == RowScope)
758 : {
759 0 : Int k=0;
760 0 : for (Int i=0;i<(pars_.doAC_ ? pars_.nAnt_ : pars_.nAnt_-1);++i) {
761 0 : for (Int j=(pars_.doAC_ ? i : i+1);j<pars_.nAnt_;++j) {
762 0 : ant2[k]=j;
763 0 : ++k;
764 : }
765 : }
766 : }
767 : else
768 0 : ant2 = thisAntenna2_;
769 0 : }
770 :
771 0 : void SimpleSimVi2::corrType (Vector<Int> & corrs) const {
772 0 : const Int& nCor(pars_.nCorr_);
773 0 : corrs.resize(nCor);
774 0 : for (Int icor=0;icor<nCor;++icor)
775 0 : corrs[icor]=corrdef_[icor];
776 0 : }
777 :
778 0 : Int SimpleSimVi2::dataDescriptionId () const { return thisSpw_; }
779 0 : void SimpleSimVi2::dataDescriptionIds (Vector<Int> & ddis) const { ddis.resize(nRows()); ddis.set(thisSpw_); }
780 0 : void SimpleSimVi2::exposure (Vector<Double> & expo) const { expo.resize(nRows()); expo.set(pars_.dt_); }
781 0 : void SimpleSimVi2::feed1 (Vector<Int> & fd1) const { fd1.resize(nRows()); fd1.set(0); }
782 0 : void SimpleSimVi2::feed2 (Vector<Int> & fd2) const { fd2.resize(nRows()); fd2.set(0); }
783 0 : void SimpleSimVi2::fieldIds (Vector<Int>& fieldids) const { fieldids.resize(nRows()); fieldids.set(thisField_); }
784 0 : void SimpleSimVi2::arrayIds (Vector<Int>& arrayids) const { arrayids.resize(nRows()); arrayids.set(0); }
785 0 : String SimpleSimVi2::fieldName () const {return "Field"+String(thisField_); }
786 :
787 0 : void SimpleSimVi2::flag (Cube<Bool> & flags) const {
788 : // unflagged
789 0 : Vector<Cube<Bool>> flagsVec;
790 0 : this->flag(flagsVec);
791 0 : flags = flagsVec[0];
792 0 : }
793 :
794 0 : void SimpleSimVi2::flag (Vector<Cube<Bool>> & flags) const {
795 : // unflagged
796 0 : flags.resize(nShapes());
797 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
798 : {
799 0 : flags[ishape].resize(pars_.nCorr_,nChannPerShape_[ishape], nRowsPerShape_[ishape]);
800 0 : flags[ishape].set(false);
801 : }
802 0 : }
803 :
804 0 : void SimpleSimVi2::flagRow (Vector<Bool> & rowflags) const { rowflags.resize(nRows()); rowflags.set(false); }
805 0 : void SimpleSimVi2::observationId (Vector<Int> & obsids) const { obsids.resize(nRows()); obsids.set(0); }
806 0 : Int SimpleSimVi2::polarizationId () const { return 0; }
807 0 : void SimpleSimVi2::processorId (Vector<Int> & procids) const { procids.resize(nRows()); procids.set(0); }
808 0 : void SimpleSimVi2::scan (Vector<Int> & scans) const { scans.resize(nRows()); scans.set(thisScan_); }
809 0 : String SimpleSimVi2::sourceName () const { return "Source"+String(thisField_); }
810 0 : void SimpleSimVi2::stateId (Vector<Int> & stateids) const { stateids.resize(nRows()); stateids.set(0); }
811 :
812 0 : Int SimpleSimVi2::polFrame () const { return 0; } // SSVi2NotYetImplemented() }
813 :
814 0 : void SimpleSimVi2::spectralWindows (Vector<Int> & spws) const
815 : {
816 0 : spws.resize(nRows());
817 0 : if(pars_.spwScope_ == ChunkScope || pars_.spwScope_ == SubchunkScope)
818 0 : spws.set(thisSpw_);
819 : else
820 : {
821 0 : indgen(spws);
822 0 : if(pars_.antennaScope_ == RowScope)
823 0 : spws = spws / nBsln_;
824 : }
825 0 : }
826 :
827 0 : void SimpleSimVi2::polarizationIds (Vector<Int> & polIds) const { polIds.resize(nRows()); polIds.set(0);}
828 0 : void SimpleSimVi2::time (Vector<Double> & t) const { t.resize(nRows()); t.set(thisTime_); }
829 0 : void SimpleSimVi2::timeCentroid (Vector<Double> & t) const { t.resize(nRows()); t.set(thisTime_); }
830 :
831 0 : void SimpleSimVi2::timeInterval (Vector<Double> & ti) const { ti.resize(nRows()); ti.set(pars_.dt_); }
832 0 : void SimpleSimVi2::uvw (Matrix<Double> & uvwmat) const { uvwmat.resize(3,nRows()); uvwmat.set(0); } // zero for now
833 :
834 0 : void SimpleSimVi2::visibilityCorrected (Cube<Complex> & vis) const {
835 : // from DATA, for now
836 0 : this->visibilityObserved(vis);
837 0 : }
838 :
839 0 : void SimpleSimVi2::visibilityCorrected (Vector<Cube<Complex>> & vis) const {
840 : // from DATA, for now
841 0 : this->visibilityObserved(vis);
842 0 : }
843 :
844 0 : void SimpleSimVi2::visibilityModel (Cube<Complex> & vis) const {
845 0 : Vector<Cube<Complex>> visVec;
846 0 : this->visibilityModel(visVec);
847 0 : vis = visVec[0];
848 0 : }
849 :
850 0 : void SimpleSimVi2::visibilityModel (Vector<Cube<Complex>> & vis) const {
851 0 : vis.resize(nShapes());
852 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
853 : {
854 0 : vis[ishape].resize(pars_.nCorr_,nChannPerShape_[ishape], nRowsPerShape_[ishape]);
855 0 : for (int icor=0;icor<pars_.nCorr_;++icor)
856 : {
857 0 : vis[ishape](Slice(icor),Slice(),Slice()).set(vis0_(icor,thisField_));
858 : }
859 : }
860 0 : }
861 :
862 0 : void SimpleSimVi2::visibilityObserved (Cube<Complex> & vis) const {
863 0 : Vector<Cube<Complex>> visVec;
864 0 : this->visibilityObserved(visVec);
865 0 : vis = visVec[0];
866 0 : }
867 :
868 0 : void SimpleSimVi2::visibilityObserved (Vector<Cube<Complex>> & vis) const {
869 0 : vis.resize(nShapes());
870 0 : this->visibilityModel(vis);
871 :
872 0 : Vector<Int> a1;
873 0 : Vector<Int> a2;
874 0 : this->antenna1(a1);
875 0 : this->antenna2(a2);
876 :
877 0 : Array<Complex> specvis;
878 0 : Matrix<Float> G(pars_.gain_);
879 0 : Matrix<Float> Tsys(pars_.tsys_);
880 :
881 0 : rownr_t vbRowIdx = 0;
882 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
883 : {
884 0 : rownr_t thisShapeVbRowIdx = vbRowIdx;
885 0 : if (pars_.doParang_ && pars_.nCorr_==4)
886 0 : corruptByParang(vis[ishape], thisShapeVbRowIdx);
887 :
888 0 : if (abs(pars_.c0_)>0.0) { // Global offset for systematic solve testing
889 0 : Cube<Complex> v(vis[ishape](Slice(0,1,1),Slice(),Slice()));
890 0 : v*=(pars_.c0_);
891 0 : if (pars_.nCorr_>1) {
892 0 : v.reference(vis[ishape](Slice(pars_.nCorr_-1,1,1),Slice(),Slice()));
893 0 : v*=(conj(pars_.c0_*pars_.c0_)); // twice the phase in the opposite direction
894 : }
895 0 : }
896 :
897 0 : for (rownr_t irow=0;irow<nRowsPerShape_[ishape];++irow, ++vbRowIdx) {
898 0 : for (int icorr=0;icorr<pars_.nCorr_;++icorr) {
899 0 : specvis.reference(vis[ishape](Slice(icorr),Slice(),Slice(irow)));
900 0 : specvis*=sqrt( G(icorr/2,a1(vbRowIdx)) * G(icorr%2,a2(vbRowIdx)) );
901 0 : if (pars_.doNorm_)
902 0 : specvis/=sqrt( Tsys(icorr/2,a1(vbRowIdx)) * Tsys(icorr%2,a2(vbRowIdx)) );
903 : }
904 : }
905 :
906 : // Now add noise
907 0 : if (pars_.doNoise_)
908 0 : this->addNoise(vis[ishape], thisShapeVbRowIdx);
909 : }
910 0 : }
911 :
912 0 : void SimpleSimVi2::floatData (Cube<Float> & fcube) const {
913 0 : Vector<Cube<Float>> fCubeVec;
914 0 : this->floatData(fCubeVec);
915 0 : fcube = fCubeVec[0];
916 0 : }
917 :
918 0 : void SimpleSimVi2::floatData (Vector<Cube<Float>> & fcubes) const {
919 0 : fcubes.resize(nShapes());
920 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
921 : {
922 0 : fcubes[ishape].resize(pars_.nCorr_,nChannPerShape_[ishape], nRowsPerShape_[ishape]);
923 : // set according to stokes
924 : // TBD
925 0 : fcubes[ishape] = 0.0;
926 : // add noise
927 : // TBD
928 : }
929 0 : }
930 :
931 0 : IPosition SimpleSimVi2::visibilityShape () const { return IPosition(3,pars_.nCorr_,pars_.nChan_(thisSpw_),nRows()); }
932 :
933 0 : void SimpleSimVi2::sigma (Matrix<Float> & sigmat) const {
934 0 : Vector<Matrix<Float>> sigVec;
935 0 : this->sigma(sigVec);
936 0 : sigmat = sigVec[0];
937 0 : }
938 :
939 0 : void SimpleSimVi2::sigma (Vector<Matrix<Float>> & sigmat) const {
940 0 : Vector<Matrix<Float>> wtMatVec;
941 0 : this->weight(wtMatVec);
942 0 : sigmat.resize(nShapes());
943 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
944 : {
945 0 : sigmat[ishape].resize(pars_.nCorr_, nRowsPerShape_[ishape]);
946 0 : sigmat[ishape] = (1.f/sqrt(wtMatVec[ishape]));
947 : }
948 0 : }
949 :
950 0 : void SimpleSimVi2::weight (Matrix<Float> & wtmat) const {
951 0 : Vector<Matrix<Float>> wtVec;
952 0 : this->weight(wtVec);
953 0 : wtmat = wtVec[0];
954 0 : }
955 :
956 0 : void SimpleSimVi2::weight (Vector<Matrix<Float>> & wtMatVec) const {
957 0 : wtMatVec.resize(nShapes());
958 :
959 0 : Vector<Int> spws;
960 0 : Vector<Int> a1;
961 0 : Vector<Int> a2;
962 0 : this->spectralWindows(spws);
963 0 : this->antenna1(a1);
964 0 : this->antenna2(a2);
965 :
966 0 : rownr_t vbRowIdx = 0;
967 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
968 : {
969 0 : wtMatVec[ishape].resize(pars_.nCorr_, nRowsPerShape_[ishape]);
970 0 : for (rownr_t irow=0;irow<nRowsPerShape_[ishape];++irow, ++vbRowIdx) {
971 0 : wtMatVec[ishape](Slice(),Slice(irow)) = wt0_(spws(vbRowIdx)); // spw-specific
972 : // non-ACs have twice this weight
973 0 : if(a1[vbRowIdx] != a2[vbRowIdx])
974 : {
975 0 : Array<Float> thiswtmat(wtMatVec[ishape](Slice(),Slice(irow)));
976 0 : thiswtmat*=Float(2.0);
977 0 : }
978 : }
979 : }
980 0 : }
981 :
982 0 : Bool SimpleSimVi2::weightSpectrumExists () const { return true; }
983 0 : Bool SimpleSimVi2::sigmaSpectrumExists () const { return true; }
984 :
985 0 : void SimpleSimVi2::weightSpectrum (Cube<Float> & wtsp) const {
986 0 : Vector<Cube<Float>> wtspVec;
987 0 : this->weightSpectrum(wtspVec);
988 0 : wtsp = wtspVec[0];
989 0 : }
990 :
991 0 : void SimpleSimVi2::weightSpectrum (Vector<Cube<Float>> & wtsp) const {
992 0 : wtsp.resize(nShapes());
993 :
994 0 : Vector<Int> spws;
995 0 : Vector<Int> a1;
996 0 : Vector<Int> a2;
997 0 : this->spectralWindows(spws);
998 0 : this->antenna1(a1);
999 0 : this->antenna2(a2);
1000 :
1001 0 : Matrix<Float> Tsys(pars_.tsys_);
1002 :
1003 0 : rownr_t vbRowIdx = 0;
1004 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
1005 : {
1006 0 : wtsp[ishape].resize(pars_.nCorr_,nChannPerShape_[ishape], nRowsPerShape_[ishape]);
1007 :
1008 0 : for (rownr_t irow=0;irow<nRowsPerShape_[ishape];++irow, ++vbRowIdx) {
1009 0 : wtsp[ishape](Slice(),Slice(),Slice(irow)) = wt0_(spws(vbRowIdx)); // spw-specific
1010 0 : if (pars_.doNoise_)
1011 : {
1012 : // non-ACs have twice this weight
1013 0 : if(a1[vbRowIdx] != a2[vbRowIdx])
1014 : {
1015 0 : Array<Float> thiswtsp(wtsp[ishape](Slice(),Slice(),Slice(irow)));
1016 0 : thiswtsp*=Float(2.0);
1017 0 : }
1018 0 : if (!pars_.doNorm_) {
1019 0 : for (Int icorr=0;icorr<pars_.nCorr_;++icorr) {
1020 0 : Array<Float> thiswt(wtsp[ishape](Slice(icorr),Slice(),Slice(irow)));
1021 0 : Float c= Tsys(icorr/2,a1(vbRowIdx))*Tsys(icorr%2,a2(vbRowIdx));
1022 0 : thiswt/=c;
1023 0 : }
1024 : }
1025 : }
1026 : }
1027 : }
1028 0 : }
1029 :
1030 0 : void SimpleSimVi2::sigmaSpectrum (Cube<Float> & sigsp) const {
1031 0 : Vector<Cube<Float>> sigspVec;
1032 0 : this->sigmaSpectrum(sigspVec);
1033 0 : sigsp = sigspVec[0];
1034 0 : }
1035 :
1036 0 : void SimpleSimVi2::sigmaSpectrum (Vector<Cube<Float>> & sigsp) const {
1037 0 : Vector<Cube<Float>> wtVec;
1038 0 : this->weightSpectrum(wtVec);
1039 0 : sigsp.resize(nShapes());
1040 0 : for (rownr_t ishape = 0 ; ishape < nShapes(); ++ishape)
1041 : {
1042 0 : sigsp[ishape].resize(pars_.nCorr_,nChannPerShape_[ishape], nRowsPerShape_[ishape]);
1043 0 : sigsp[ishape] = (1.f/sqrt(wtVec[ishape]));
1044 : }
1045 0 : }
1046 :
1047 0 : const casacore::Vector<casacore::Float>& SimpleSimVi2::feed_pa (casacore::Double t) const
1048 : {
1049 0 : feedpa_.resize(nAntennas(),0.0f);
1050 0 : if (pars_.doParang_)
1051 0 : feedpa_.set(0.05f*Float(t-t0_)); // 0.05 rad/s
1052 :
1053 0 : return feedpa_;
1054 : }
1055 :
1056 :
1057 : // +=========================+
1058 : // | |
1059 : // | Chunk and MS Level Data |
1060 : // | |
1061 : // +=========================+
1062 :
1063 :
1064 0 : CountedPtr<WeightScaling> SimpleSimVi2::getWeightScaling () const {
1065 0 : return WeightScaling::generateUnityWeightScaling();
1066 : }
1067 :
1068 0 : MEpoch SimpleSimVi2::getEpoch () const { SSVi2NotYetImplemented() }
1069 :
1070 0 : Vector<Int> SimpleSimVi2::getCorrelations () const {
1071 : // The correlation indices
1072 0 : Vector<Int> corrs(pars_.nCorr_);
1073 0 : indgen(corrs);
1074 0 : return corrs;
1075 0 : }
1076 :
1077 0 : Vector<Int> SimpleSimVi2::getChannels (Double, Int, Int spw, Int) const {
1078 0 : Vector<Int> chans(pars_.nChan_(spw));
1079 0 : indgen(chans);
1080 0 : return chans;
1081 0 : }
1082 :
1083 0 : Vector<Double> SimpleSimVi2::getFrequencies (Double, Int, Int spw, Int) const {
1084 0 : return pars_.freqs(spw);
1085 : /*
1086 : Vector<Double> freqs(pars_.nChan_(spw));
1087 : indgen(freqs);
1088 : freqs*=pars_.df_(spw);
1089 : freqs+=(pars_.df_(spw)/2. + pars_.refFreq_(spw));
1090 : return freqs;
1091 : */
1092 : }
1093 :
1094 0 : Vector<Double> SimpleSimVi2::getChanWidths (Double, Int, Int spw, Int) const {
1095 0 : Vector<Double> widths(pars_.nChan_(spw),pars_.df_(spw));
1096 0 : return widths;
1097 : }
1098 :
1099 : // get back the selected spectral windows and spectral channels for
1100 : // current ms
1101 :
1102 0 : const SpectralWindowChannels & SimpleSimVi2::getSpectralWindowChannels (Int /*msId*/, Int /*spectralWindowId*/) const { SSVi2NotYetImplemented() }
1103 :
1104 : // Return number of antennasm spws, polids, ddids
1105 :
1106 0 : Int SimpleSimVi2::nAntennas () const { return pars_.nAnt_; }
1107 0 : Int SimpleSimVi2::nDataDescriptionIds () const { return pars_.nSpw_; }
1108 0 : Int SimpleSimVi2::nPolarizationIds () const { return Int(1); }
1109 :
1110 0 : rownr_t SimpleSimVi2::nRowsInChunk () const { SSVi2NotYetImplemented() } // number rows in current chunk
1111 0 : rownr_t SimpleSimVi2::nRowsViWillSweep () const { SSVi2NotYetImplemented() } // number of rows in all selected ms's
1112 :
1113 0 : Int SimpleSimVi2::nSpectralWindows () const { return pars_.nSpw_; }
1114 :
1115 0 : Int SimpleSimVi2::nTimes() const {
1116 0 : SSVi2NotYetImplemented();
1117 : }
1118 :
1119 :
1120 0 : void SimpleSimVi2::configureNewSubchunk() {
1121 0 : nShapes_ = 1;
1122 0 : rownr_t spwRowFactor = 1, antennaRowFactor = 1;
1123 :
1124 0 : if(pars_.spwScope_ == RowScope)
1125 0 : spwRowFactor = pars_.nSpw_;
1126 0 : if(pars_.antennaScope_ == RowScope)
1127 0 : antennaRowFactor = nBsln_;
1128 :
1129 0 : if(pars_.spwScope_ == RowScope)
1130 0 : nShapes_ = pars_.nSpw_;
1131 :
1132 0 : nRows_ = spwRowFactor * antennaRowFactor;
1133 :
1134 0 : nRowsPerShape_.resize(nShapes_);
1135 0 : nRowsPerShape_ = nRows_ / nShapes_;
1136 :
1137 0 : Vector<Int> spws;
1138 0 : spectralWindows(spws);
1139 0 : nChannPerShape_.resize(nShapes_);
1140 0 : if(pars_.spwScope_ == RowScope)
1141 0 : nChannPerShape_ = pars_.nChan_;
1142 : else
1143 0 : nChannPerShape_ = pars_.nChan_(spws(0));
1144 0 : nCorrsPerShape_.resize(nShapes_);
1145 0 : nCorrsPerShape_ = pars_.nCorr_;
1146 :
1147 : // Poke the vb to do this
1148 0 : vb_->configureNewSubchunk(0,"faked",false,
1149 0 : isNewArrayId(),isNewFieldId(),
1150 0 : isNewSpectralWindow(),getSubchunkId(),
1151 0 : nRowsPerShape_,
1152 0 : nChannPerShape_,
1153 0 : nCorrsPerShape_,
1154 0 : getCorrelations(),
1155 0 : corrdef_,corrdef_,
1156 0 : WeightScaling::generateUnityWeightScaling());
1157 0 : }
1158 :
1159 : // Generate noise on data
1160 0 : void SimpleSimVi2::addNoise(Cube<Complex>& vis, rownr_t vbRowOffset) const {
1161 :
1162 0 : IPosition sh3(vis.shape());
1163 :
1164 0 : Int64 seed(thisScan_*1000000+thisField_*100000+thisSpw_*10000 + Int(thisTime_-t0_));
1165 0 : ACG r(seed);
1166 0 : Normal rn(&r,0.0,1.0/wt0_(thisSpw_));
1167 :
1168 0 : Vector<Int> a1;
1169 0 : Vector<Int> a2;
1170 0 : this->antenna1(a1);
1171 0 : this->antenna2(a2);
1172 :
1173 0 : Matrix<Float> Tsys(pars_.tsys_);
1174 0 : for (Int i=0;i<sh3[2];++i) {
1175 0 : Float c(1.0);
1176 0 : if (a1(i+vbRowOffset)!=a2(i+vbRowOffset))
1177 0 : c=1./sqrt(2.0);
1178 0 : for (Int k=0;k<sh3[0];++k) {
1179 0 : Float d(c);
1180 0 : if (!pars_.doNorm_)
1181 0 : d*=sqrt(Tsys(k/2,a1(i+vbRowOffset))*Tsys(k%2,a2(i+vbRowOffset)));
1182 0 : for (Int j=0;j<sh3[1];++j) {
1183 0 : vis(k,j,i)+=Complex(d)*Complex(rn(),rn());
1184 : }
1185 : }
1186 : }
1187 :
1188 0 : }
1189 :
1190 0 : void SimpleSimVi2::corruptByParang(Cube<Complex>& vis, rownr_t vbRowOffset) const {
1191 :
1192 :
1193 : // Assumes constant time in this subchunk
1194 0 : Vector<Float> pa(this->feed_pa(thisTime_));
1195 : // Row-wise antenna Ids
1196 0 : Vector<Int> a1;
1197 0 : Vector<Int> a2;
1198 0 : this->antenna1(a1);
1199 0 : this->antenna2(a2);
1200 :
1201 0 : Cube<Complex> v;
1202 0 : if (pars_.polBasis_=="circ") {
1203 0 : Matrix<Complex> P(2,pars_.nAnt_,Complex(1.0f));
1204 0 : for (Int iant=0;iant<pars_.nAnt_;++iant) {
1205 0 : Float& a(pa(iant));
1206 0 : P(1,iant)=Complex(cos(a),sin(a));
1207 0 : P(0,iant)=conj(P(1,iant));
1208 : }
1209 0 : for (Int irow=0;irow<vis.shape()(2);++irow) {
1210 0 : Int& i(a1(irow+vbRowOffset)), j(a2(irow+vbRowOffset));
1211 0 : for (Int icor=0;icor<pars_.nCorr_;++icor) {
1212 0 : v.reference(vis(Slice(icor,1,1),Slice(),Slice(irow,1,1)));
1213 0 : v*=(P(icor/2,i)*conj(P(icor%2,j))); // cf PJones
1214 : }
1215 : }
1216 0 : }
1217 0 : else if (pars_.polBasis_=="lin") {
1218 0 : Matrix<Float> P(2,pars_.nAnt_,0.0f);
1219 0 : for (Int iant=0;iant<pars_.nAnt_;++iant) {
1220 0 : Float& a(pa(iant));
1221 0 : P(0,iant)=cos(a);
1222 0 : P(1,iant)=sin(a);
1223 : }
1224 0 : Vector<Complex> v0(pars_.nCorr_),v1;
1225 0 : for (Int irow=0;irow<vis.shape()(2);++irow) {
1226 0 : Int& i(a1(irow+vbRowOffset)), j(a2(irow+vbRowOffset));
1227 0 : Float& C(P(0,i)),S(P(1,i)),c(P(0,j)),s(P(1,j));
1228 0 : for (Int ich=0;ich<pars_.nChan_(thisSpw_);++ich) {
1229 0 : v1.reference(vis.xyPlane(irow).column(ich));
1230 0 : v0.assign(v1); // _deep copy_ of this irow
1231 0 : v1(0)=( C*v0[0]+S*v0[2])*c + ( C*v0[1]+S*v0[3])*s;
1232 0 : v1(1)=( C*v0[0]+S*v0[2])*(-s) + ( C*v0[1]+S*v0[3])*c;
1233 0 : v1(2)=(-S*v0[0]+C*v0[2])*c + (-S*v0[1]+C*v0[3])*s;
1234 0 : v1(3)=(-S*v0[0]+C*v0[2])*(-s) + (-S*v0[1]+C*v0[3])*c;
1235 : }
1236 : }
1237 0 : }
1238 0 : }
1239 :
1240 0 : const casacore::MSAntennaColumns& SimpleSimVi2::antennaSubtablecols() const
1241 : {
1242 0 : return *antennaSubTablecols_p;
1243 : }
1244 :
1245 : // Access to dataDescription subtable
1246 0 : const casacore::MSDataDescColumns& SimpleSimVi2::dataDescriptionSubtablecols() const
1247 : {
1248 0 : return *ddSubTablecols_p;
1249 : }
1250 :
1251 : // Access to feed subtable
1252 0 : const casacore::MSFeedColumns& SimpleSimVi2::feedSubtablecols() const
1253 : {
1254 0 : SSVi2NotYetImplemented();
1255 : }
1256 :
1257 : // Access to field subtable
1258 0 : const casacore::MSFieldColumns& SimpleSimVi2::fieldSubtablecols() const
1259 : {
1260 0 : SSVi2NotYetImplemented();
1261 : }
1262 :
1263 : // Access to flagCmd subtable
1264 0 : const casacore::MSFlagCmdColumns& SimpleSimVi2::flagCmdSubtablecols() const
1265 : {
1266 0 : SSVi2NotYetImplemented();
1267 : }
1268 :
1269 : // Access to history subtable
1270 0 : const casacore::MSHistoryColumns& SimpleSimVi2::historySubtablecols() const
1271 : {
1272 0 : SSVi2NotYetImplemented();
1273 : }
1274 :
1275 : // Access to observation subtable
1276 0 : const casacore::MSObservationColumns& SimpleSimVi2::observationSubtablecols() const
1277 : {
1278 0 : SSVi2NotYetImplemented();
1279 : }
1280 :
1281 : // Access to pointing subtable
1282 0 : const casacore::MSPointingColumns& SimpleSimVi2::pointingSubtablecols() const
1283 : {
1284 0 : SSVi2NotYetImplemented();
1285 : }
1286 :
1287 : // Access to polarization subtable
1288 0 : const casacore::MSPolarizationColumns& SimpleSimVi2::polarizationSubtablecols() const
1289 : {
1290 0 : return *polSubTablecols_p;
1291 : }
1292 :
1293 : // Access to processor subtable
1294 0 : const casacore::MSProcessorColumns& SimpleSimVi2::processorSubtablecols() const
1295 : {
1296 0 : SSVi2NotYetImplemented();
1297 : }
1298 :
1299 : // Access to spectralWindow subtable
1300 0 : const casacore::MSSpWindowColumns& SimpleSimVi2::spectralWindowSubtablecols() const
1301 : {
1302 0 : return *spwSubTablecols_p;
1303 : }
1304 :
1305 : // Access to state subtable
1306 0 : const casacore::MSStateColumns& SimpleSimVi2::stateSubtablecols() const
1307 : {
1308 0 : SSVi2NotYetImplemented();
1309 : }
1310 :
1311 : // Access to doppler subtable
1312 0 : const casacore::MSDopplerColumns& SimpleSimVi2::dopplerSubtablecols() const
1313 : {
1314 0 : SSVi2NotYetImplemented();
1315 : }
1316 :
1317 : // Access to freqOffset subtable
1318 0 : const casacore::MSFreqOffsetColumns& SimpleSimVi2::freqOffsetSubtablecols() const
1319 : {
1320 0 : SSVi2NotYetImplemented();
1321 : }
1322 :
1323 : // Access to source subtable
1324 0 : const casacore::MSSourceColumns& SimpleSimVi2::sourceSubtablecols() const
1325 : {
1326 0 : SSVi2NotYetImplemented();
1327 : }
1328 :
1329 : // Access to sysCal subtable
1330 0 : const casacore::MSSysCalColumns& SimpleSimVi2::sysCalSubtablecols() const
1331 : {
1332 0 : SSVi2NotYetImplemented();
1333 : }
1334 :
1335 : // Access to weather subtable
1336 0 : const casacore::MSWeatherColumns& SimpleSimVi2::weatherSubtablecols() const
1337 : {
1338 0 : SSVi2NotYetImplemented();
1339 : }
1340 :
1341 0 : SimpleSimVi2Factory::SimpleSimVi2Factory(const SimpleSimVi2Parameters& pars)
1342 0 : : pars_(pars)
1343 0 : {}
1344 :
1345 0 : SimpleSimVi2Factory::~SimpleSimVi2Factory () {}
1346 :
1347 0 : ViImplementation2 * SimpleSimVi2Factory::createVi () const {
1348 :
1349 0 : ViImplementation2* vii = new SimpleSimVi2(pars_);
1350 0 : return vii;
1351 :
1352 : }
1353 :
1354 0 : SimpleSimVi2LayerFactory::SimpleSimVi2LayerFactory(const SimpleSimVi2Parameters& pars)
1355 : : ViiLayerFactory(),
1356 0 : pars_(pars)
1357 : {
1358 : // cout << "&pars = " << &pars << endl;
1359 : // cout << "&pars_ = " << &pars_ << endl;
1360 0 : }
1361 :
1362 :
1363 : // SimpleSimVi2-specific layer-creater
1364 0 : ViImplementation2 * SimpleSimVi2LayerFactory::createInstance (ViImplementation2* /*vii0*/) const {
1365 :
1366 : // No deeper layers
1367 : // Assert(vii0==NULL);
1368 :
1369 : // Make it and return it
1370 0 : ViImplementation2 *vii = new SimpleSimVi2(pars_);
1371 0 : return vii;
1372 : }
1373 :
1374 :
1375 : } // end namespace vi
1376 :
1377 : } //# NAMESPACE CASA - END
1378 :
1379 :
|
