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