Line data Source code
1 : // Based on code/alma/apps/UVWCoords
2 :
3 : #include <casacore/casa/Logging/LogIO.h>
4 : #include <msvis/MSVis/MSUVWGenerator.h>
5 : #include <casacore/measures/Measures/MEpoch.h>
6 : #include <casacore/measures/Measures/MFrequency.h>
7 : #include <casacore/measures/Measures/MPosition.h>
8 : #include <casacore/ms/MeasurementSets/MSColumns.h>
9 : #include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
10 : #include <casacore/measures/Measures/MCBaseline.h>
11 : #include <casacore/casa/Utilities/GenSort.h>
12 : #include <map>
13 :
14 : using namespace casacore;
15 : namespace casa {
16 :
17 : // The UvwCoords ctor has lines for the antennas, antenna offsets, and station
18 : // positions. This ctor assumes they're present in msc_p if present at all.
19 0 : MSUVWGenerator::MSUVWGenerator(MSColumns &msc_ref, const MBaseline::Types bltype,
20 0 : const Muvw::Types) :
21 : msc_p(msc_ref),
22 : bl_csys_p(MBaseline::Ref(bltype)), // MBaseline::J2000, ITRF, etc.
23 0 : antColumns_p(msc_p.antenna()),
24 0 : antPositions_p(antColumns_p.positionMeas()),
25 0 : antOffset_p(antColumns_p.offsetMeas()),
26 0 : refpos_p(antPositions_p(0)), // We use the first antenna for the reference
27 0 : feedOffset_p(msc_p.feed().positionMeas())
28 : {
29 : // It seems that using a String is the only safe way to go from say,
30 : // MPosition::ITRF to MBaseline::ITRF.
31 0 : MBaseline::getType(refposref_p,
32 0 : MPosition::showType(refpos_p.getRef().getType()));
33 :
34 0 : fill_bl_an(bl_an_p);
35 0 : }
36 :
37 0 : MSUVWGenerator::~MSUVWGenerator(){
38 0 : }
39 :
40 0 : void MSUVWGenerator::fill_bl_an(Vector<MVBaseline>& bl_an_p)
41 : {
42 0 : nant_p = antPositions_p.table().nrow();
43 0 : logSink() << LogIO::DEBUG1 << "nant_p: " << nant_p << LogIO::POST;
44 :
45 0 : Double max_baseline = -1.0;
46 : Double bl_len;
47 :
48 0 : const ScalarColumn<Double>& antDiams = antColumns_p.dishDiameter();
49 0 : Double smallestDiam = antDiams(0);
50 0 : Double secondSmallestDiam = antDiams(0);
51 :
52 0 : bl_an_p.resize(nant_p);
53 0 : for(uInt an = 0; an < nant_p; ++an){
54 : // MVBaselines are basically xyz Vectors, not Measures.
55 0 : bl_an_p[an] = MVBaseline(refpos_p.getValue(), antPositions_p(an).getValue());
56 :
57 : // MVBaseline has functions to return the length, but Manhattan distances
58 : // are good enough for this, and faster than a sqrt.
59 0 : Vector<Double> bluvw(bl_an_p[an].getValue());
60 0 : bl_len = fabs(bluvw[0]) + fabs(bluvw[1]) + fabs(bluvw[2]);
61 :
62 0 : if(bl_len > max_baseline)
63 0 : max_baseline = bl_len;
64 :
65 0 : if(antDiams(an) < secondSmallestDiam){
66 0 : if(antDiams(an) < smallestDiam){
67 0 : secondSmallestDiam = smallestDiam;
68 0 : smallestDiam = antDiams(an);
69 : }
70 : else
71 0 : secondSmallestDiam = antDiams(an);
72 : }
73 : }
74 :
75 : // Setting timeRes_p to 0.0025 * the time for a 1 radian phase change on the
76 : // longest baseline at 2x the primary beamwidth should be sufficiently short
77 : // for Earth based observations. Space-based baselines will move faster, but
78 : // probably don't have the data rate to support full beam imaging. An
79 : // alternative limit could come from |d(uvw)/dt|/|uvw|, but that is
80 : // guaranteed to be at least somewhat larger than this limit (replace the
81 : // 2.44 with 2, and remove the dish diameters and max_baseline).
82 : //
83 : // Do not raise the 0.0025 coefficient by too much, since the times used for
84 : // UVW calculation could be biased by as much as -0.5 * timeRes_p if more
85 : // than one integration fits into timeRes_p. timeRes_p is not intended so
86 : // much for skipping calculations as it is for preventing antUVW_p from being
87 : // calculated for each row in the same integration. The integration interval
88 : // may change within an MS, but ideally antUVW_p is calculated once per
89 : // integration (timeRes_p <~ integration) and there is no bias.
90 0 : timeRes_p = 0.0025 * 24.0 * 3600.0 / (6.283 * 2.44) *
91 0 : sqrt(smallestDiam * secondSmallestDiam) / max_baseline;
92 0 : }
93 :
94 0 : void MSUVWGenerator::uvw_an(const MEpoch& timeCentroid, const Int fldID,
95 : const Bool WSRTConvention)
96 : {
97 0 : const MDirection& phasedir = msc_p.field().phaseDirMeas(fldID);
98 0 : MeasFrame measFrame(refpos_p, timeCentroid, phasedir);
99 0 : MVBaseline mvbl;
100 0 : MBaseline basMeas;
101 :
102 : //logSink() << LogIO::DEBUG1
103 : // << "timeCentroid: " << timeCentroid
104 : // << "\nfldID: " << fldID
105 : // << "\nphasedir: " << phasedir
106 : // << LogIO::POST;
107 :
108 0 : MBaseline::Ref basref(refposref_p);
109 0 : basMeas.set(mvbl, basref); // in antenna frame
110 0 : basMeas.getRefPtr()->set(measFrame);
111 :
112 : // Setup a machine for converting a baseline vector from the antenna frame to
113 : // bl_csys_p's frame
114 0 : MBaseline::Convert elconv(basMeas, bl_csys_p);
115 :
116 0 : for(uInt i = 0; i < nant_p; ++i){
117 : //TODO: (Soon!) Antenna offsets are not handled yet.
118 0 : basMeas.set(bl_an_p[i], basref);
119 0 : MBaseline basOutFrame = elconv(basMeas);
120 : //MBaseline::Types botype = MBaseline::castType(basOutFrame.getRef().getType());
121 0 : MVuvw uvwOutFrame(basOutFrame.getValue(), phasedir.getValue());
122 :
123 0 : antUVW_p[i] = uvwOutFrame.getValue();
124 : }
125 :
126 : // Tony Willis supplied this crucial bit of history: the WSRT definition of
127 : // phase is opposite to that of the VLA - the WSRT definition of the l,m
128 : // direction cosines is that l increases toward decreasing RA, whereas the
129 : // VLA definition is that l increases toward increasing RA. AIPS seems to
130 : // understand this. (For a completely useless piece of trivia, the WSRT defn
131 : // of phase is consistent with that of the Cambridge one-mile telescope as
132 : // Wim Brouw set things up to be in agreement with the one-mile telescope ca
133 : // 1968.)
134 : // RR: It is not just l, v and w are flipped as well.
135 0 : if(WSRTConvention)
136 0 : antUVW_p = -antUVW_p;
137 0 : }
138 :
139 : // antUVW_p must be set up for the correct timeCentroid and phase direction by
140 : // uvw_an() before calling this.
141 0 : void MSUVWGenerator::uvw_bl(const uInt ant1, const uInt,// feed1,
142 : const uInt ant2, const uInt,// feed2,
143 : Array<Double>& uvw)
144 : {
145 : //uvw.resize(3); // Probably redundant. Does it significantly slow things down?
146 : //TODO: Feed offsets are not handled yet.
147 0 : uvw = antUVW_p[ant2] - antUVW_p[ant1];
148 0 : }
149 :
150 0 : Bool MSUVWGenerator::make_uvws(const Vector<Int> flds)
151 : {
152 0 : ArrayColumn<Double>& UVWcol = msc_p.uvw();
153 0 : const ScalarMeasColumn<MEpoch>& timeCentMeas = msc_p.timeCentroidMeas();
154 0 : const ScalarColumn<Int> fieldID(msc_p.fieldId());
155 0 : const ScalarColumn<Int> ant1(msc_p.antenna1());
156 0 : const ScalarColumn<Int> ant2(msc_p.antenna2());
157 0 : const ScalarColumn<Int> feed1(msc_p.feed1());
158 0 : const ScalarColumn<Int> feed2(msc_p.feed2());
159 0 : const ScalarColumn<Int> obsID(msc_p.observationId());
160 :
161 : // Use a time ordered index to minimize the number of calls to uvw_an.
162 : // TODO: use field as a secondary sort key.
163 0 : Vector<uInt> tOI;
164 0 : GenSortIndirect<Double>::sort(tOI, msc_p.timeCentroid().getColumn());
165 :
166 : // Having uvw_an() calculate positions for each antenna for every field is
167 : // somewhat inefficient since in a multiconfig MS not all antennas will be
168 : // used in each time interval, but it's not clear that determining which
169 : // antennas will be used for a given uvw_an() call would be any more
170 : // efficient. It's not horribly inefficient, because uvw_an() is O(nant_p),
171 : // and uvw_bl() is only called for baselines that are actually used.
172 0 : antUVW_p.resize(nant_p);
173 :
174 0 : logSink() << LogOrigin("MSUVWGenerator", "make_uvws") << LogIO::NORMAL3;
175 :
176 0 : logSink() << LogIO::DEBUG1 << "timeRes_p: " << timeRes_p << LogIO::POST;
177 :
178 : try{
179 0 : Bool oldWsrtConvention = (msc_p.observation().telescopeName()(obsID(tOI[0])) ==
180 : "WSRT");
181 :
182 : // Ensure a call to uvw_an on the 1st iteration.
183 0 : const Unit sec("s");
184 0 : Double oldTime = timeCentMeas(tOI[0]).get(sec).getValue() - 2.0 * timeRes_p;
185 0 : Int oldFld = -2;
186 :
187 : // IO reordering buffer to avoid poor IO patterns on files not sorted by time
188 0 : std::map<uInt, Vector<Double> > writebuffer;
189 :
190 0 : for(uInt row = 0; row < msc_p.nrow(); ++row){
191 0 : uInt toir = tOI[row];
192 0 : const MEpoch & timecentmeas = timeCentMeas(toir);
193 0 : Double currTime = timecentmeas.get(sec).getValue();
194 0 : Int currFld = fieldID(toir);
195 0 : Bool newWsrtConvention = (msc_p.observation().telescopeName()(obsID(toir)) ==
196 : "WSRT");
197 0 : if ((row % (1<<22)) == 0) {
198 : logSink() << LogIO::NORMAL << "Rows handled: "
199 0 : << row << "/" << msc_p.nrow() << LogIO::POST;
200 : }
201 :
202 0 : if(currTime - oldTime > timeRes_p || currFld != oldFld
203 0 : || newWsrtConvention != oldWsrtConvention){
204 0 : oldTime = currTime;
205 0 : oldFld = currFld;
206 :
207 0 : if(newWsrtConvention != oldWsrtConvention){
208 : logSink() << LogIO::WARN
209 : << "There is a switch to or from the WSRT convention at row "
210 : << toir << ".\nWatch for an inconsistency in the sign of UVW."
211 0 : << LogIO::POST;
212 0 : oldWsrtConvention = newWsrtConvention;
213 : }
214 :
215 : //logSink() << LogIO::DEBUG1 << "currTime: " << currTime
216 : // << "\ncurrFld: " << currFld << LogIO::POST;
217 0 : uvw_an(timecentmeas, currFld, newWsrtConvention);
218 : }
219 :
220 0 : if(flds[fieldID(toir)] > -1){
221 : // uvw_bl(ant1(toir), ant2(toir),
222 : // feed1(toir), feed2(toir), UVWcol(toir));
223 : // collect output and write larger chunk in row sequential order
224 0 : writebuffer[toir] = antUVW_p[ant2(toir)] - antUVW_p[ant1(toir)];
225 0 : if (writebuffer.size() > 100000) {
226 0 : for (auto it = writebuffer.begin(); it != writebuffer.end(); ++it) {
227 0 : UVWcol.put(it->first, it->second);
228 : }
229 0 : writebuffer.clear();
230 : }
231 : }
232 : }
233 : // flush rest of buffer
234 0 : for (auto it = writebuffer.begin(); it != writebuffer.end(); ++it) {
235 0 : UVWcol.put(it->first, it->second);
236 : }
237 : }
238 0 : catch(AipsError x){
239 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
240 0 : << LogIO::POST;
241 0 : throw(AipsError("Error in MSUVWGenerator::make_uvws."));
242 : return false;
243 : }
244 0 : return true;
245 : }
246 :
247 : // void MSUVWGenerator::get_ant_offsets(const MDirection& dir_with_a_frame)
248 : // {
249 : // // This appears to be a required column of the ANTENNA table in version 2.0
250 : // // of the MeasurementSet definition
251 : // // (http://aips2.nrao.edu/docs/notes/229/229.html), so it is assumed to be
252 : // // present. However, it is usually a set of zeroes, based on the common
253 : // // belief that it is only needed for heterogeneous arrays, since the
254 : // // receivers of homogeneous arrays move in concert. That is not true when
255 : // // there are independent pointing errors.
256 :
257 : // // Convert ant_offset_measures to Vectors and check for nonzeroness.
258 : // ant_offset_vec_p.resize(nant);
259 : // for(uInt n = 0; n < nant; ++n)
260 : // ant_offset_vec_p[n] = ant_offset_meas_p.convert(0, pointingdir).getValue();
261 : // }
262 :
263 : // Bool MSUVWGenerator::set_receiv_offsets(const MDirection& dir_with_a_frame)
264 : // {
265 : // // This appears to be a required column of the FEED table in version 2.0
266 : // // of the MeasurementSet definition
267 : // // (http://aips2.nrao.edu/docs/notes/229/229.html), so it is assumed to be
268 : // // present. However, it is usually a set of zeroes, based on the common
269 : // // belief that it is only needed for heterogeneous arrays, since the
270 : // // receivers of homogeneous arrays move in concert. That is not true when
271 : // // there are independent pointing errors.
272 :
273 : // // Convert feed_offset_measures to Vectors and check for nonzeroness.
274 : // Vector<Vector<Double> > offsetvects;
275 : // offsetvects.resize(nant);
276 : // for(uInt n = 0; n < nant; ++n){
277 : // offsetvects[n] = ant_offsets->convert(0, pointingdir).getValue();
278 : // if(ant_offsets[n] != ant_offsets[0]){
279 : // varying_offsets = true;
280 : // }
281 : // }
282 :
283 : // ignore_offsets = true;
284 : // return ignore_offsets;
285 : // }
286 :
287 : using namespace casacore;
288 : } // Ends namespace casa.
|
