Line data Source code
1 : //# BeamSkyJones.h: Definitions of interface for BeamSkyJones
2 : //# Copyright (C) 1996,1997,1998,1999,2000,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 adressed 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 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_TRANSFORM2_BEAMSKYJONES_H
30 : #define SYNTHESIS_TRANSFORM2_BEAMSKYJONES_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <casacore/casa/Containers/Block.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 : #include <casacore/measures/Measures/MDirection.h>
36 : #include <casacore/measures/Measures.h>
37 : #include <casacore/measures/Measures/Stokes.h>
38 : #include <synthesis/TransformMachines2/SkyJones.h>
39 : #include <synthesis/TransformMachines/PBMath.h>
40 : #include <msvis/MSVis/VisBufferUtil.h>
41 : #include <casacore/casa/Arrays/ArrayFwd.h>
42 :
43 : namespace casacore{
44 :
45 : class ImageRegion;
46 : class String;
47 :
48 : }
49 :
50 : namespace casa {
51 : //#forward
52 : //# Need forward declaration for Solve in the Jones Matrices
53 :
54 : class SkyModel;
55 :
56 : namespace refim { //# namespace for refactoring
57 :
58 :
59 : // <summary> beam-like sky-plane effects for the SkyEquation </summary>
60 :
61 : // <use visibility=export>
62 :
63 : // <reviewed reviewer="" date="" tests="" demos="">
64 :
65 : // <prerequisite>
66 : // <li> <linkto class="SkyEquation">SkyEquation</linkto> class
67 : // <li> <linkto class="SkyJones">SkyJones</linkto> class
68 : // <li> <linkto class="PBMathInterface">PBMathInterface</linkto> class
69 : // </prerequisite>
70 : //
71 : // <etymology>
72 : // BeamSkyJones, derived from SkyJones, describes an interface to
73 : // beam-based SkyJones objects. Like SkyJones, it too is an Abstract
74 : // Base Class, but implements the beam-related methods.
75 : // </etymology>
76 : //
77 : // <synopsis>
78 : //
79 : // <motivation>
80 : // The properties of sky-plane based calibration effects must be described
81 : // for the <linkto class="SkyEquation">SkyEquation</linkto>; this class
82 : // encapsulates the antenna beam-based aspects which are present in at least
83 : // a few other specific SkyJones (VPSkyJones and DBeamSkyJones).
84 : // </motivation>
85 : //
86 : // <todo asof="98/11/01">
87 : // <li> Solvable part needs implementation: we need to derive an
88 : // image of gradients of the elements of the Jones matrix. See VisJones
89 : // for how to do this.
90 : // <li> The casacore::MS, version II, will have a beam subtable, which will have PBMaths
91 : // for each antenna. Until this becomes available, we need to do some
92 : // fudging; for example, we
93 : // </todo>
94 :
95 : class BeamSkyJones : virtual public SkyJones {
96 :
97 : public:
98 :
99 : // Eventually, the casacore::MS will have all the beam information in its Beam Subtable.
100 : // Till then, we either guess the PB to use or explicitly define it upon construction.
101 : // Construct from a Measurement Set, figure out the most appropriate PBMath object
102 : // from casacore::MS information
103 : BeamSkyJones(const casacore::Quantity ¶llacticAngleIncrement = casacore::Quantity(720.0, "deg"), // def= 1 PA interval
104 : BeamSquint::SquintType doSquint = BeamSquint::NONE, // def= no beam squint offsets
105 : const casacore::Quantity &skyPositionThreshold = casacore::Quantity(180,"deg")); // def= assume there is no change of
106 : // this operator due to position offset
107 :
108 :
109 : // Virtual destructor (this is a virtual base class)
110 : virtual ~BeamSkyJones() = 0;
111 :
112 : // Print out information concerning the state of this object
113 : virtual void showState(casacore::LogIO& os);
114 :
115 : // Apply Jones matrix to an image (and adjoint)
116 : // No "applyInverse" is available from the SkyJones classes,
117 : // you can get them directly from PBMath or you can get
118 : // the equivalent effect by dividing by grad grad Chi^2
119 : // in ImageSkyModel.
120 : // <group>
121 : casacore::ImageInterface<casacore::Complex>& apply(const casacore::ImageInterface<casacore::Complex>& in,
122 : casacore::ImageInterface<casacore::Complex>& out,
123 : const vi::VisBuffer2& vb, casacore::Int row,
124 : casacore::Bool forward=true);
125 : casacore::ImageInterface<casacore::Float>& apply(const casacore::ImageInterface<casacore::Float>& in,
126 : casacore::ImageInterface<casacore::Float>& out,
127 : const vi::VisBuffer2& vb, casacore::Int row);
128 :
129 : casacore::ImageInterface<casacore::Float>& applySquare(const casacore::ImageInterface<casacore::Float>& in,
130 : casacore::ImageInterface<casacore::Float>& out,
131 : const vi::VisBuffer2& vb, casacore::Int row);
132 : // </group>
133 :
134 : // Apply Jones matrix to a sky component (and adjoint)
135 : // <group>
136 : SkyComponent& apply(SkyComponent& in,
137 : SkyComponent& out,
138 : const vi::VisBuffer2& vb, casacore::Int row,
139 : casacore::Bool forward = true);
140 : SkyComponent& applySquare(SkyComponent& in,
141 : SkyComponent& out,
142 : const vi::VisBuffer2& vb, casacore::Int row);
143 : // </group>
144 :
145 : // Understand if things have changed since last PB application
146 : // <group>
147 :
148 : // Reset to pristine state
149 : virtual void reset();
150 :
151 : // Has this operator changed since the last Application?
152 : // (or more properly, since the last update() )
153 : virtual casacore::Bool changed(const vi::VisBuffer2& vb, casacore::Int row);
154 :
155 : // Does the operator change in this visbuffer or since the last call?
156 : // May not be useful -- check it out: m.a.h. Dec 30 1999
157 : virtual casacore::Bool change(const vi::VisBuffer2& vb);
158 :
159 : // Does this operator changed in this VisBuffer,
160 : // starting with row1?
161 : // If yes, we return in row2, the last row that has the
162 : // same SkyJones as row1.
163 : // NOTE: need to first call changed(const VisBuffer& vb, casacore::Int row) and
164 : // shield the user from the case where the fieldID has changed
165 : // (which only changes in blocks)
166 : virtual casacore::Bool changedBuffer(const vi::VisBuffer2& vb, casacore::Int row1, casacore::Int& row2);
167 :
168 : // Update the FieldID, Telescope, pointingDirection, Parallactic angle info
169 : void update(const vi::VisBuffer2& vb, casacore::Int row);
170 :
171 : // if (changed) {update()}
172 : virtual void assure (const vi::VisBuffer2& vb, casacore::Int row);
173 : // </group>
174 :
175 : // Return the type of this Jones matrix (actual type of derived class).
176 0 : virtual casa::SkyJones::Type type() {return casa::SkyJones::E;};
177 :
178 : // Apply gradient
179 : virtual casacore::ImageInterface<casacore::Complex>&
180 : applyGradient(casacore::ImageInterface<casacore::Complex>& result, const vi::VisBuffer2& vb,
181 : casacore::Int row);
182 :
183 : virtual SkyComponent&
184 : applyGradient(SkyComponent& result, const vi::VisBuffer2& vb,
185 : casacore::Int row);
186 :
187 : // Is this solveable?
188 0 : virtual casacore::Bool isSolveable() {return false;};
189 :
190 : // Initialize for gradient search
191 : virtual void initializeGradients();
192 :
193 : // Finalize for gradient search
194 : virtual void finalizeGradients();
195 :
196 : // Add to Gradient Chisq
197 : virtual void addGradients(const vi::VisBuffer2& vb, casacore::Int row, const casacore::Float sumwt,
198 : const casacore::Float chisq, const casacore::Matrix<casacore::Complex>& c,
199 : const casacore::Matrix<casacore::Float>& f);
200 :
201 : // Solve
202 : //virtual casacore::Bool solve (SkyEquation& se);
203 :
204 : // Manage the PBMath objects
205 : // <group>
206 : // set the PB based on telescope name, antennaID and feedID
207 : // If antennaID or feedID is -1, the PBMath object is set for
208 : // all antennae or feeds, respectively. These are the default
209 : // values to retain the previous interface.
210 : //
211 : // Note. It would be nice to change the interface and make antennaID
212 : // and feedID the second and the third parameter, respectively.
213 : void setPBMath(const casacore::String &telescope, PBMath &myPBmath,
214 230 : const casacore::Int &antennaID = -1, const casacore::Int &feedID = -1);
215 :
216 : // get the PBMath object; returns false if that one doesn't exist,
217 : // true if it does exist and is OK
218 : // whichAnt is an index into an array of PBMaths, which is different
219 : // for all telescope/antenna/feed combinations
220 : // Not sure why we need such a low-level method declared as public,
221 : // retained to preserve old interface
222 : casacore::Bool getPBMath(casacore::uInt whichAnt, PBMath &myPBMath) const;
223 :
224 : // get the PBMath object; returns false if that one doesn't exist,
225 : // true if it does exist and is OK
226 : // antennaID and feedID default to -1 to preserve the old interface
227 : // TODO: change the interface and make antennaID and feedID the
228 : // second and third parameter respectively to have a better looking code
229 : casacore::Bool getPBMath(const casacore::String &telescope, PBMath &myPBMath,
230 : const casacore::Int &antennaID = -1, const casacore::Int &feedID = -1) const;
231 :
232 : casacore::Quantity getPAIncrement() {return casacore::Quantity(parallacticAngleIncrement_p,"rad");}
233 :
234 : casacore::Quantity getSkyPositionThreshold() {return casacore::Quantity(skyPositionThreshold_p,"rad");}
235 :
236 : // Return true if all antennas share a common VP
237 : casacore::Bool isHomogeneous() const;
238 : //</group>
239 :
240 : // Get the casacore::ImageRegion of the primary beam on an Image for a given pointing
241 : // Note: casacore::ImageRegion is not necesarily constrained to lie within the
242 : // image region (for example, if the pointing center is near the edge of the
243 : // image). fPad: extra padding over the primary beam supporrt,
244 : // fractional (ie, 1.2 for 20% padding), in all directions.
245 : // (note: we do not properly treat squint yet, this will cover it for now)
246 : // iChan: frequency channel to take: lowest frequency channel is safe for all
247 : //
248 : // Potential problem: this casacore::ImageRegion includes all casacore::Stokes and Frequency Channels
249 : // present in the input image.
250 : //COMMENTING out for now as this depend on PBMathInterface and which depends
251 : //back on SkyJones::sizeType
252 : casacore::ImageRegion* extent (const casacore::ImageInterface<casacore::Complex>& im,
253 : const vi::VisBuffer2& vb,
254 : const casacore::Int irow=-1,
255 : const casacore::Float fPad=1.2,
256 : const casacore::Int iChan=0,
257 : const casa::SkyJones::SizeType sizeType=casa::SkyJones::COMPOSITE);
258 :
259 : casacore::ImageRegion* extent (const casacore::ImageInterface<casacore::Float>& im,
260 : const vi::VisBuffer2& vb, const casacore::Int irow=-1,
261 : const casacore::Float fPad=1.2, const casacore::Int iChan=0,
262 : const casa::SkyJones::SizeType sizeType=casa::SkyJones::COMPOSITE);
263 :
264 : // summarize the PBMaths contained here.
265 : // n = -1 => terse table
266 : // n = 0 => table plus constructor values
267 : // n = m => plot m samples of the PB profile
268 : virtual void summary(casacore::Int n=0);
269 :
270 : //return the telescope it is on at this state
271 : casacore::String telescope();
272 :
273 : //Get an idea of the support of the PB in number of pixels
274 : virtual casacore::Int support(const vi::VisBuffer2& vb, const casacore::CoordinateSystem& cs);
275 :
276 : private:
277 :
278 :
279 : casacore::String telescope_p;
280 :
281 : casacore::Int lastFieldId_p;
282 :
283 : casacore::Int lastArrayId_p;
284 :
285 : casacore::Int lastMSId_p;
286 : casacore::Double lastTime_p;
287 : BeamSquint::SquintType doSquint_p;
288 :
289 : casacore::Double parallacticAngleIncrement_p; // a parallactic angle threshold
290 : // beyond which the operator is considered to be
291 : // changed (in radians)
292 : casacore::Double skyPositionThreshold_p; // a sky position threshold beyond
293 : // which the operator is considered to be changed
294 : // (in radians)
295 : casacore::Block<casacore::Double> lastParallacticAngles_p; // a cache of parallactic angles
296 : // used when the operator was applied last time.
297 : // One value in radians for each beam model in PBMaths.
298 : // A zero-length block means that the operator
299 : // has never been applied
300 : casacore::Block<casacore::MDirection> lastDirections_p; // a chache of directions
301 : // used when the operator was applied last time.
302 : // One element for each beam model in PBMaths.
303 : // A zero-length block means that the operator
304 : // has never been applied
305 :
306 : // One or more PBMaths (a common one for the
307 : // entire array, or one for each antenna)
308 : // This requires some sorting out for heterogeneous arrays!
309 : casacore::Block<PBMath> myPBMaths_p;
310 : // Names of telescopes (parralel with PBMaths
311 : casacore::Block<casacore::String> myTelescopes_p;
312 :
313 : // Antenna IDs (parallel with PBMaths)
314 : casacore::Block<casacore::Int> myAntennaIDs_p;
315 : // Feed IDs (parallel with PBMaths)
316 : casacore::Block<casacore::Int> myFeedIDs_p;
317 :
318 : // cache of the indices to the PBMaths container for antenna/feed 1 and 2
319 : mutable casacore::CountedPtr<vi::VisBuffer2> lastUpdateVisBuffer_p; // to ensure that the cache
320 : // is filled for the correct
321 : // VisBuffer. The value is used
322 : // for comparison only
323 : mutable casacore::Int lastUpdateRow_p; // to ensure that the cache is filled for
324 : // correct row in the VisBuffer
325 : mutable casacore::Int lastUpdateIndex1_p; // index of the first antenna/feed
326 : mutable casacore::Int lastUpdateIndex2_p; // index of the second antenna/feed
327 : //
328 :
329 : mutable casacore::Bool hasBeenApplied; // true if the operator has been applied at least once
330 :
331 : // update the indices cache. This method could be made protected in the
332 : // future (need access functions for lastUpdateIndex?_p to benefit from it)
333 : // Cache will be valid for a given VisBuffer and row
334 : void updatePBMathIndices(const vi::VisBuffer2 &vb, casacore::Int row) const;
335 :
336 : protected:
337 : // return true if two directions are close enough to consider the
338 : // operator unchanged, false otherwise
339 : casacore::Bool directionsCloseEnough(const casacore::MDirection &dir1,
340 : const casacore::MDirection &dir2) const;
341 :
342 : // return index of compareTelescope, compareAntenna and compareFeed in
343 : // myTelescopes_p, myAntennaIDs and myFeedIDs; -1 if not found
344 : // if compareAntenna or compareTelescope is -1, this means that the
345 : // PBMath class is the same for all antennae/feeds. An exception will
346 : // be raised, if separate PBMath objects have been assigned by setPBMath
347 : // for different feeds/antennae but -1 is used for query.
348 : //
349 : // It would be good to rename this function to indexBeams as this name
350 : // is more appropriate.
351 : //
352 : casacore::Int indexTelescope(const casacore::String & compareTelescope,
353 0 : const casacore::Int &compareAntenna=-1,
354 213 : const casacore::Int &compareFeed=-1) const;
355 : casacore::MDirection convertDir(const vi::VisBuffer2& vb, const casacore::MDirection& inDir, const casacore::MDirection::Types outType);
356 : casacore::CountedPtr<VisBufferUtil> vbutil_p;
357 :
358 : };
359 :
360 : } //# end of namespace refim
361 :
362 : } // end namespace casa
363 :
364 : #endif
365 :
|