Line data Source code
1 : //# SimACohCalc.cc: Simulated additive errors
2 : //# Copyright (C) 1996,1997,1999,2000,2001
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 adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #include <synthesis/MeasurementComponents/SimACohCalc.h>
30 : #include <msvis/MSVis/VisBuffer.h>
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
33 : #include <casacore/casa/Logging/LogIO.h>
34 : #include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
35 : #include <casacore/measures/Measures.h>
36 :
37 : using namespace casacore;
38 : namespace casa { //# NAMESPACE CASA - BEGIN
39 :
40 : // Note: this simplistic implementation just generates a new random
41 : // noise value for every call of apply, it doesn't keep track of
42 : // time and antenna to return the same value if called with the same
43 : // coordinates. Thus applyInv will only work correctly if
44 : // called from a separate object in the exact same sequence as apply.
45 :
46 0 : SimACohCalc::SimACohCalc(const Int seed,
47 : const Float antefficiency,
48 : const Float correfficiency,
49 : const Float spillefficiency,
50 : const Float tau,
51 : const Quantity& trx,
52 : const Quantity& tatmos,
53 0 : const Quantity& tcmb) :
54 : rndGen_p(seed),
55 : // noiseDist_p(&rndGen_p, 0.0, 0.5),
56 : noiseDist_p(&rndGen_p, 0.0, 1.0),
57 : antefficiency_p(antefficiency),
58 : correfficiency_p(correfficiency),
59 : spillefficiency_p(spillefficiency),
60 : tau_p(tau),
61 : trx_p(trx),
62 : tatmos_p(tatmos),
63 0 : tcmb_p(tcmb)
64 : {
65 0 : }
66 :
67 :
68 0 : SimACohCalc::~SimACohCalc()
69 0 : {};
70 :
71 :
72 0 : VisBuffer& SimACohCalc::apply(VisBuffer& vb)
73 : {
74 0 : Complex c[4];
75 :
76 : // Double averageNoise = 0.0;
77 : // Int nAverage = 0;
78 : // Double averageT = 0.0;
79 : // Int nAverage = 0;
80 :
81 : // In case you are confused, the 1e-4 converts the Diam from cm to m
82 0 : Double fact = 4 * C::sqrt2 * 1.38062e-16 * 1e23 * 1e-4 /
83 0 : ( antefficiency_p * correfficiency_p * C::pi );
84 :
85 0 : Double lasttime = 0.0;
86 0 : Double time = 0.0;
87 0 : Vector<MDirection> azel;
88 : Int ant1, ant2;
89 0 : Double el1=0.0, el2=0.0, airmass1=0.0, airmass2=0.0;
90 0 : Double tsys=0.0, sigma = 1.0;
91 : Double diam1, diam2;
92 : Float deltaNu;
93 : Double tint;
94 0 : Vector<Double> antDiams;
95 :
96 0 : Double trx_k = trx_p.getValue("K");
97 0 : Double tatmos_k = tatmos_p.getValue("K");
98 0 : Double tcmb_k = tcmb_p.getValue("K");
99 0 : Bool zeroSpacing = false;
100 :
101 0 : LogIO os(LogOrigin("SimACohCalc", "apply()", WHERE));
102 :
103 : {
104 0 : Int iSpW=vb.spectralWindow();
105 0 : deltaNu = vb.msColumns().spectralWindow().totalBandwidth()(iSpW) /
106 0 : Float(vb.msColumns().spectralWindow().numChan()(iSpW));
107 :
108 : // live dangerously: assume all vis have the same tint
109 0 : tint = vb.msColumns().exposure()(0);
110 :
111 0 : antDiams = vb.msColumns().antenna().dishDiameter().getColumn();
112 : }
113 :
114 0 : for (Int row=0; row<vb.nRow(); row++) {
115 :
116 0 : ant1=vb.antenna1()(row);
117 0 : ant2=vb.antenna2()(row);
118 0 : diam1 = antDiams(ant1);
119 0 : diam2 = antDiams(ant2);
120 :
121 0 : time = vb.time()(row);
122 0 : if (time != lasttime) {
123 0 : azel = vb.azel(time);
124 : }
125 0 : el1 = azel(ant1).getAngle("rad").getValue()(1);
126 0 : el2 = azel(ant2).getAngle("rad").getValue()(1);
127 :
128 0 : airmass1 = 1.0; airmass2 = 1.0;
129 0 : if ( (el1 > 0.0 && el2 > 0.0) || tau_p == 0.0) {
130 0 : airmass1 = 1.0/ sin(el1);
131 0 : airmass2 = 1.0/ sin(el2);
132 :
133 0 : tsys = trx_k * sqrt(exp(tau_p * airmass1)) * sqrt(exp(tau_p * airmass2))
134 0 : + tatmos_k *
135 0 : (sqrt(exp(tau_p * airmass1))*sqrt(exp(tau_p * airmass2)) - spillefficiency_p)
136 : + tcmb_k;
137 :
138 0 : sigma = fact * tsys / diam1 / diam2 / sqrt( deltaNu * tint );
139 :
140 0 : vb.sigma()(row) = sigma;
141 0 : sigma = max( sigma, 1e-9 );
142 0 : vb.weight()(row) *= 1.0/ pow(sigma, 2.0);
143 :
144 :
145 0 : zeroSpacing = false;
146 0 : if (vb.uvw()(row)(0) == 0.0 && vb.uvw()(row)(1) == 0.0) {
147 0 : zeroSpacing = true;
148 : }
149 :
150 0 : for (Int chn=0; chn<vb.nChannel(); chn++) {
151 0 : for (Int i=0; i<4; i++) {
152 0 : if (zeroSpacing) {
153 0 : Float re = 1.41421356*noiseDist_p();
154 0 : c[i]= Float(sigma) * Complex(re, 0.0);
155 : } else {
156 : // huh. why would one use a Bessel function (product of normals) here?
157 : // Float re = noiseDist_p()*noiseDist_p();
158 : // c[i]= Float(sigma) * Complex(re);
159 0 : c[i]= Float(sigma) * Complex(noiseDist_p(),noiseDist_p());
160 : }
161 : // nAverage++; averageNoise += sigma;
162 : }
163 :
164 0 : CStokesVector noiseCoh(c);
165 0 : vb.visibility()(chn,row)+=noiseCoh;
166 : }
167 : // nAverage++; averageT += tsys;
168 : }
169 : }
170 : // if (nAverage > 0) {
171 : // averageNoise /= Float(nAverage);
172 : // os << "Average noise added to visibilities: " << averageNoise << " Jy" << LogIO::POST;
173 : // }
174 : // if (nAverage > 0) {
175 : // averageT /= Float(nAverage);
176 : // os << "Noise added to visibilities with average Tsys=" << averageT << " K" << LogIO::POST;
177 : // }
178 0 : return vb;
179 : };
180 :
181 :
182 0 : VisBuffer& SimACohCalc::applyInv(VisBuffer& vb)
183 : {
184 0 : Complex c[4];
185 0 : for (Int row=0; row<vb.nRow(); row++) {
186 0 : for (Int chn=0; chn<vb.nChannel(); chn++) {
187 0 : for (Int i=0; i<4; i++) c[i]=Complex(noiseDist_p(),noiseDist_p());
188 0 : CStokesVector noiseCoh(c);
189 0 : vb.visibility()(chn,row)-=noiseCoh;
190 : }
191 : }
192 0 : return vb;
193 : }
194 :
195 : } //# NAMESPACE CASA - END
196 :
|
