Line data Source code
1 : //# BeamSkyJones.cc: Implementation for BeamSkyJones
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 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 <casacore/casa/aips.h>
30 : #include <casacore/casa/BasicSL/Complex.h>
31 :
32 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
33 : #include <casacore/ms/MeasurementSets/MSColumns.h>
34 : #include <casacore/ms/MeasurementSets/MSObsColumns.h>
35 : #include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
36 : #include <casacore/tables/Tables.h>
37 : #include <casacore/measures/Measures/Stokes.h>
38 : #include <casacore/measures/Measures/MeasConvert.h>
39 :
40 : #include <casacore/casa/BasicSL/Constants.h>
41 : #include <casacore/measures/Measures/MeasTable.h>
42 : #include <components/ComponentModels/Flux.h>
43 : #include <components/ComponentModels/ComponentShape.h>
44 :
45 : #include <synthesis/TransformMachines2/BeamSkyJones.h>
46 : #include <synthesis/TransformMachines/PBMath.h>
47 :
48 : #include <msvis/MSVis/VisBuffer2.h>
49 : #include <msvis/MSVis/VisibilityIterator2.h>
50 : #include <msvis/MSVis/VisBufferComponents2.h>
51 : #include <casacore/images/Images/ImageInterface.h>
52 : #include <casacore/images/Regions/ImageRegion.h>
53 :
54 : #include <casacore/casa/Utilities/Assert.h>
55 :
56 :
57 : using namespace casacore;
58 : namespace casa{
59 :
60 :
61 : namespace refim{
62 : // temporary, for debugging
63 : /*#include <casa/Quanta/MVAngle.h>
64 : void printDirection(std::ostream &os,const MDirection &dir) throw (AipsError) {
65 : double lngbuf=dir.getValue().getLong("deg").getValue();
66 : if (lngbuf<0) lngbuf+=360.;
67 : os<<(dir.getRefString()!="GALACTIC"?MVAngle::Format(MVAngle::TIME):
68 : MVAngle::Format(MVAngle::ANGLE))<<MVAngle(Quantity(lngbuf,"deg"))<<" "
69 : <<MVAngle(dir.getValue().getLat("deg"))<<
70 : " ("<<dir.getRefString()<<")";
71 : }
72 : */
73 : //
74 :
75 : using namespace casa;
76 : using namespace casa::refim;
77 :
78 215 : BeamSkyJones::BeamSkyJones(
79 : const Quantity ¶llacticAngleIncrement,
80 : BeamSquint::SquintType doSquint,
81 0 : const Quantity &skyPositionThreshold) :
82 : doSquint_p(doSquint),
83 215 : parallacticAngleIncrement_p(parallacticAngleIncrement.getValue("rad")),
84 215 : skyPositionThreshold_p(skyPositionThreshold.getValue("rad")),
85 : lastUpdateVisBuffer_p( ), lastUpdateRow_p(-1),
86 645 : lastUpdateIndex1_p(-1), lastUpdateIndex2_p(-1), hasBeenApplied(false), vbutil_p(nullptr)
87 :
88 : {
89 215 : reset();
90 215 : setThreshold(0.01); // use this in apply to determine level of cutoff
91 215 : };
92 :
93 215 : void BeamSkyJones::reset()
94 : {
95 215 : lastFieldId_p=-1;
96 215 : lastArrayId_p=-1;
97 215 : lastMSId_p=0;
98 215 : lastTime_p=-1.0;
99 215 : telescope_p="";
100 215 : }
101 :
102 197 : BeamSkyJones::~BeamSkyJones()
103 : {
104 197 : };
105 :
106 0 : void BeamSkyJones::showState(LogIO& os)
107 : {
108 0 : os << "Field ID = " << lastFieldId_p+1 << LogIO::POST;
109 0 : os << "Telescope = " << telescope_p << LogIO::POST;
110 0 : for (uInt i=0;i<lastParallacticAngles_p.nelements();++i) {
111 : os << "ParAngle[d] ("<<i<<" model) = " <<
112 0 : lastParallacticAngles_p[i]/C::pi*180.<< LogIO::POST;
113 : os << "Pointing direction ("<<i<<" model) = "<<
114 0 : lastDirections_p[i].getAngle().getValue("deg")(0) <<
115 0 : ", " << lastDirections_p[i].getAngle().getValue("deg")(1) << LogIO::POST;
116 : }
117 0 : os << "delta PA[d] = " << Quantity(parallacticAngleIncrement_p,"rad").getValue("deg") << LogIO::POST;
118 0 : os << "skyPositionThreshold[d] = " << Quantity(skyPositionThreshold_p,"rad").getValue("deg") << LogIO::POST;
119 0 : os << "SquintType = " << (Int)(doSquint_p) << LogIO::POST;
120 0 : };
121 :
122 4 : String BeamSkyJones::telescope(){
123 4 : return telescope_p;
124 : }
125 : // update the indices cache. This method could be made protected in the
126 : // future (need access functions for lastUpdateIndex?_p to benefit from it)
127 : // Cache will be valid for a given VisBuffer and row
128 96 : void BeamSkyJones::updatePBMathIndices(const vi::VisBuffer2 &vb, Int row) const
129 : {
130 : // for debug
131 : // cout<<endl<<"BeamSkyJones::updatePBMathIndices row="<<row<<endl<<endl;
132 : //
133 :
134 : // we will not check msId, arrayId and fieldID as they are supposed to be
135 : // checked before this method is called. However, if this method is to
136 : // be made protected, a change may be required here to avoid nasty
137 : // surprises in the future.
138 :
139 :
140 : ///////////////AS visbuffer2 copy does not seem to work and there is no way to know the state of the previous vb...will update it anyways.
141 : /*if (!lastUpdateVisBuffer_p.null() && vb.rowIds()(0)==lastUpdateVisBuffer_p->rowIds()(0) && row==lastUpdateRow_p) return;
142 : lastUpdateVisBuffer_p=vi::VisBuffer2::factory(vi::VbRekeyable);
143 : // lastUpdateVisBuffer_p->setShape(vb.getShape()(0), vb.getShape()(1), vb.getShape()(2));
144 : //lastUpdateVisBuffer_p->copy(vb,true);
145 :
146 : vi::VisBufferComponents2 comps = vi::VisBufferComponents2::these(vi::VisBufferComponent2::DataDescriptionIds, vi::VisBufferComponent2::FieldId, vi::VisBufferComponent2::RowIds);
147 : lastUpdateVisBuffer_p->copyComponents(vb, comps, true, true);
148 : */
149 96 : lastUpdateRow_p=row;
150 :
151 : // Getting the reference on antennae/feed IDs is a
152 : // fast operation as caching is implemented inside VisBuffer.
153 96 : DebugAssert(row<(Int)vb.antenna1().nelements(),AipsError);
154 96 : Int ant1=vb.antenna1()[row];
155 96 : Int ant2=vb.antenna2()[row];
156 96 : Int feed1=vb.feed1()[row];
157 96 : Int feed2=vb.feed2()[row];
158 :
159 : // telescope_p should be valid at this stage because it is updated
160 : // after each ArrayID change. Care must be taken if the method is to be
161 : // made protected.
162 96 : lastUpdateIndex1_p=indexTelescope(telescope_p,ant1,feed1);
163 96 : lastUpdateIndex2_p=indexTelescope(telescope_p,ant2,feed2);
164 96 : }
165 :
166 96 : Bool BeamSkyJones::changed(const vi::VisBuffer2& vb, Int row)
167 : {
168 : // for debug
169 : // cout<<endl<<"BeamSkyJones::changed row="<<row<<" lastUpdateRow_p="<<
170 : // lastUpdateRow_p<<endl<<endl;
171 : //
172 :
173 96 : if (row < 0) row = 0;
174 :
175 154 : if(vb.msId() != lastMSId_p || vb.arrayId()(0)!=lastArrayId_p ||
176 58 : vb.fieldId()(0)!=lastFieldId_p) {
177 38 : lastUpdateVisBuffer_p=NULL; // invalidate index cache
178 38 : return true;
179 : }
180 :
181 58 : MDirection::Types pointingdirType=MDirection::castType(lastDirections_p[lastUpdateIndex1_p].getRef().getType());
182 : //if in a local frame and time change => pointing most probably changed
183 58 : if((pointingdirType >= MDirection::HADEC) && (pointingdirType <= MDirection::AZELSWGEO)){
184 58 : if(lastTime_p != vb.time()(0))
185 0 : return True;
186 : }
187 :
188 : //if (lastUpdateIndex1_p<0 || lastUpdateIndex2_p<0) return true;
189 :
190 58 : updatePBMathIndices(vb,row); // lastUpdateIndex?_p are now valid
191 :
192 : //Unnecessary ...i believe and causes issues with PSF making
193 : //if (!hasBeenApplied) return true; // we shouldn't have such a flag in
194 : // a well designed code
195 :
196 58 : if (!lastParallacticAngles_p.nelements() && myPBMaths_p.nelements())
197 0 : return true; // it's a first call of this method and setPBMath has
198 : // definitely been called before
199 :
200 :
201 :
202 : // Obtaining a reference on parallactic angles is a fast operation as
203 : // caching is implemented inside VisBuffer.
204 58 : Float feed1_pa=vb.feedPa1()[row];
205 58 : Float feed2_pa=vb.feedPa2()[row];
206 :
207 : // it may be good to check here whether an indexed beam model
208 : // depend on parallactic angle before returning true
209 : // An additional interface function may be required for PBMath classes
210 :
211 58 : if (lastUpdateIndex1_p!=-1)
212 58 : if (abs(feed1_pa-lastParallacticAngles_p[lastUpdateIndex1_p]) >
213 58 : parallacticAngleIncrement_p) return true;
214 :
215 58 : if (lastUpdateIndex2_p!=-1)
216 58 : if (abs(feed2_pa-lastParallacticAngles_p[lastUpdateIndex2_p]) >
217 58 : parallacticAngleIncrement_p) return true;
218 :
219 : /*
220 : These direction test are not used right now and are terrible calculations to do on
221 : billions of rows of data
222 : If it is needed a faster algorithm for changed direction is needed
223 : ////
224 : if (lastUpdateIndex1_p!=-1)
225 : if (!directionsCloseEnough(lastDirections_p[lastUpdateIndex1_p],
226 : vb.direction1()[row])) return true;
227 :
228 : if (lastUpdateIndex2_p!=-1)
229 : if (!directionsCloseEnough(lastDirections_p[lastUpdateIndex2_p],
230 : vb.direction2()[row])) return true;
231 : */
232 :
233 58 : return false;
234 : };
235 :
236 : // return true if two directions are close enough to consider the
237 : // operator unchanged, false otherwise
238 38 : Bool BeamSkyJones::directionsCloseEnough(const MDirection &dir1,
239 : const MDirection &dir2) const throw(AipsError)
240 : {
241 : Double sep;
242 38 : if (dir1.getRef()!=dir2.getRef())
243 60 : sep=dir1.getValue().separation(MDirection::Convert(dir2.getRef(),
244 60 : dir1.getRef())(dir2).getValue());
245 18 : else sep=dir1.getValue().separation(dir2.getValue());
246 38 : return (fabs(sep)<skyPositionThreshold_p);
247 : }
248 :
249 : // Does this BeamSkyJones change during this buffer, starting from
250 : // row1? If row2 <0, row2 = nRow()-1
251 0 : Bool BeamSkyJones::changedBuffer(const vi::VisBuffer2& vb, Int row1, Int& row2)
252 : {
253 0 : Int irow = row1;
254 0 : if (irow < 0) irow = 0;
255 0 : Int jrow = row2;
256 0 : if (jrow < 0) jrow = vb.nRows()-1;
257 0 : DebugAssert(jrow<vb.nRows(),AipsError);
258 :
259 : // It is not important now to have a separate function for a "block"
260 : // operation. Because an appropriate caching is implemented inside
261 : // both VisBuffer and this class, changed(vb,row) can be called in the
262 : // loop without a perfomance penalty. We preserve this function to keep
263 : // the interface intact.
264 :
265 0 : for (Int ii=irow+1;ii<=jrow;++ii)
266 0 : if (changed(vb,ii)) {
267 0 : row2 = ii-1;
268 0 : return true;
269 : }
270 0 : return false;
271 : };
272 :
273 : // as it is stated in BeamSkyJones.h this method may not be useful
274 : // anymore. Implementing it just in case it is used somewhere.
275 : // Because an appropriate caching is implemented in both VisBuffer and
276 : // BeamSkyJones, this method can use BeamSkyJones::changed in a loop
277 0 : Bool BeamSkyJones::change(const vi::VisBuffer2& vb)
278 : {
279 0 : for (rownr_t i=0;i<vb.nRows();++i)
280 0 : if (changed(vb,i)) return true;
281 0 : return false;
282 : };
283 :
284 38 : void BeamSkyJones::update(const vi::VisBuffer2& vb, Int row)
285 : {
286 : // for debug
287 : //cout<<endl<<"BeamSkyJones::update nrow="<<vb.nRow()<<" row="<<row<<" feed1="<<vb.feed1()(0)<<" feed2="<<vb.feed2()(0)<<endl<<endl;
288 : //
289 :
290 38 : if (row<0) row=0;
291 :
292 38 : lastFieldId_p=vb.fieldId()(0);
293 38 : lastArrayId_p=vb.arrayId()(0);
294 38 : lastMSId_p=vb.msId();
295 38 : lastTime_p=vb.time()(0);
296 38 : if(!vbutil_p){
297 38 : vbutil_p=new VisBufferUtil(vb);
298 : }
299 : // The pointing direction depends on feed, antenna, pointing errors, etc
300 : //MDirection pointingDirection1 = vb.direction1()(row);
301 : //MDirection pointingDirection2 = vb.direction2()(row);
302 76 : MDirection pointingDirection1 = vbutil_p->getPointingDir(vb, vb.antenna1()(row), row);
303 76 : MDirection pointingDirection2 = vbutil_p->getPointingDir(vb, vb.antenna2()(row), row);
304 : //cerr << "DIR1 " << pointingDirection1.toString() << " " << pointingDirection2.toString() << endl;
305 : // Look up correct telescope
306 38 : const MSObservationColumns& msoc=vb.subtableColumns().observation();
307 38 : telescope_p = msoc.telescopeName()(vb.arrayId()(0));
308 :
309 38 : updatePBMathIndices(vb,row); // lastUpdateIndex?_p are now valid
310 :
311 38 : if (!lastParallacticAngles_p.nelements() && myPBMaths_p.nelements()) {
312 38 : lastParallacticAngles_p.resize(myPBMaths_p.nelements());
313 38 : lastParallacticAngles_p.set(1000.); // to force recalculation
314 : // it will recalculate directions as well
315 : }
316 38 : if (!lastDirections_p.nelements() && myPBMaths_p.nelements())
317 38 : lastDirections_p.resize(myPBMaths_p.nelements());
318 :
319 76 : if (lastUpdateIndex1_p == lastUpdateIndex2_p &&
320 38 : !directionsCloseEnough(pointingDirection1,pointingDirection2)) {
321 : // the case is inhomogeneous: pointing directions are slightly
322 : // different at different antennae
323 : //This check is an overkill for standard arrays...need to find a better one
324 :
325 : /* LogIO os;
326 : os << LogIO::WARN << LogOrigin("BeamSkyJones","update")
327 : << "The pointing directions differ for different stations."
328 : << LogIO::POST << LogIO::WARN << LogOrigin("BeamSkyJones","update")
329 : << "This case is not handled correctly. Continuing anyway."<<LogIO::POST;
330 :
331 : */
332 : // we could, in principle, clone a PBMath object for one of the
333 : // antennae and rebuild lastDirections_p.
334 : // For now, the value for the second antenna will be used
335 : }
336 :
337 38 : if (lastUpdateIndex1_p!=-1)
338 38 : lastDirections_p[lastUpdateIndex1_p]=pointingDirection1;
339 : ////CAS-6688 using antenna1 only for now
340 : // if (lastUpdateIndex2_p!=-1)
341 : // lastDirections_p[lastUpdateIndex2_p]=pointingDirection2;
342 :
343 : // Obtaining a reference on parallactic angles is a fast operation as
344 : // caching is implemented inside VisBuffer.
345 38 : Float feed1_pa=vb.feedPa1()[row];
346 38 : Float feed2_pa=vb.feedPa2()[row];
347 :
348 76 : if (lastUpdateIndex1_p == lastUpdateIndex2_p &&
349 38 : abs(abs(feed1_pa-feed2_pa)-parallacticAngleIncrement_p)> 1e-5 ) {
350 : // the array is not compact: position angles are significantly
351 : // different at different antennae
352 38 : LogIO os;
353 : //Commenting out this message...more pest than useful to have it at this low level
354 : // os << LogIO::WARN << LogOrigin("BeamSkyJones","update")
355 : // << "The array is not compact, position angles differ for different stations."
356 : // << LogIO::POST << LogIO::WARN << LogOrigin("BeamSkyJones","update")
357 : // << "Primary beams are not correctly handled if they are asymmetric. Continuing anyway."<<LogIO::POST;
358 : // we could, in principle, clone a PBMath object for one of the
359 : // antennae and rebuild lastParallacticAngles_p.
360 : // For now, the value for the second antenna will be used
361 : }
362 38 : if (lastUpdateIndex1_p!=-1)
363 38 : lastParallacticAngles_p[lastUpdateIndex1_p]=feed1_pa;
364 38 : if (lastUpdateIndex2_p!=-1)
365 38 : lastParallacticAngles_p[lastUpdateIndex2_p]=feed2_pa;
366 38 : };
367 :
368 0 : void BeamSkyJones::assure (const vi::VisBuffer2& vb, Int row)
369 : {
370 0 : if(changed(vb, row)) update(vb, row);
371 0 : };
372 :
373 :
374 : ImageInterface<Complex>&
375 38 : BeamSkyJones::apply(const ImageInterface<Complex>& in,
376 : ImageInterface<Complex>& out,
377 : const vi::VisBuffer2& vb, Int row,
378 : Bool forward)
379 : {
380 38 : if(changed(vb, row)) update(vb, row);
381 38 : hasBeenApplied=true;
382 : // now lastUpdateIndex?_p are valid
383 :
384 38 : if (lastUpdateIndex1_p!=lastUpdateIndex2_p)
385 0 : throw(AipsError("BeamSkyJones::apply(Image...) - can only treat homogeneous PB case"));
386 : else {
387 : // for debug
388 : // cout<<endl<<"BeamSkyJones::apply(Image...) index="<<lastUpdateIndex1_p<<" feed="<<vb.feed1()(0)<<" direction=";
389 : // printDirection(cout,lastDirections_p[lastUpdateIndex1_p]); cout<<endl<<endl;
390 : //
391 76 : CoordinateSystem cs=in.coordinates();
392 38 : Int coordIndex=cs.findCoordinate(Coordinate::DIRECTION);
393 38 : MDirection::Types dirType=cs.directionCoordinate(coordIndex).directionType();
394 38 : PBMath myPBMath;
395 38 : if (getPBMath(lastUpdateIndex1_p, myPBMath))
396 38 : return myPBMath.applyPB(in, out, convertDir(vb, lastDirections_p[lastUpdateIndex1_p], dirType),
397 76 : Quantity(lastParallacticAngles_p[lastUpdateIndex1_p],"rad"),
398 114 : doSquint_p, false, threshold(), forward);
399 : else
400 0 : throw(AipsError("BeamSkyJones::apply(Image...)!!! - PBMath not found"));
401 : }
402 : };
403 :
404 :
405 : ImageInterface<Float>&
406 0 : BeamSkyJones::apply(const ImageInterface<Float>& in,
407 : ImageInterface<Float>& out,
408 : const vi::VisBuffer2& vb, Int row){
409 0 : if(changed(vb, row)) update(vb, row);
410 0 : hasBeenApplied=true;
411 : // now lastUpdateIndex?_p are valid
412 :
413 0 : if (lastUpdateIndex1_p!=lastUpdateIndex2_p)
414 0 : throw(AipsError("BeamSkyJones::apply(Image...) - can only treat homogeneous PB case"));
415 : else {
416 0 : PBMath myPBMath;
417 0 : CoordinateSystem cs=in.coordinates();
418 0 : Int coordIndex=cs.findCoordinate(Coordinate::DIRECTION);
419 0 : MDirection::Types dirType=cs.directionCoordinate(coordIndex).directionType();
420 0 : if (getPBMath(lastUpdateIndex1_p, myPBMath))
421 0 : return myPBMath.applyPB(in, out,convertDir(vb, lastDirections_p[lastUpdateIndex1_p], dirType),
422 0 : Quantity(lastParallacticAngles_p[lastUpdateIndex1_p],"rad"),
423 0 : doSquint_p, threshold());
424 : else
425 0 : throw(AipsError("BeamSkyJones::apply(Image...)!!! - PBMath not found"));
426 : }
427 :
428 : }
429 : ImageInterface<Float>&
430 40 : BeamSkyJones::applySquare(const ImageInterface<Float>& in,
431 : ImageInterface<Float>& out,
432 : const vi::VisBuffer2& vb, Int row)
433 : {
434 40 : if(changed(vb, row)) update(vb, row);
435 40 : hasBeenApplied=true;
436 : // now lastUpdateIndex?_p are valid
437 :
438 40 : if (lastUpdateIndex1_p!=lastUpdateIndex2_p)
439 0 : throw(AipsError("BeamSkyJones::applySquare(Image...) - can only treat homogeneous PB case"));
440 : else {
441 : // for debug
442 : //cout<<endl<<"BeamSkyJones::applySquare(Image...) index="<<lastUpdateIndex1_p<<" feed="<<vb.feed1()(0)<<" direction=";
443 : //printDirection(cout,lastDirections_p[lastUpdateIndex1_p]); cout<<endl<<endl;
444 : //
445 80 : PBMath myPBMath;
446 40 : CoordinateSystem cs=in.coordinates();
447 40 : Int coordIndex=cs.findCoordinate(Coordinate::DIRECTION);
448 40 : MDirection::Types dirType=cs.directionCoordinate(coordIndex).directionType();
449 40 : if (getPBMath(lastUpdateIndex1_p, myPBMath))
450 40 : return myPBMath.applyPB2(in, out, convertDir(vb, lastDirections_p[lastUpdateIndex1_p], dirType),
451 120 : lastParallacticAngles_p[lastUpdateIndex1_p], doSquint_p, threshold()*threshold());
452 : else
453 0 : throw(AipsError("BeamSkyJones::applySquare(Image...) - PBMath not found"));
454 : }
455 : };
456 :
457 :
458 : SkyComponent&
459 0 : BeamSkyJones::apply(SkyComponent& in,
460 : SkyComponent& out,
461 : const vi::VisBuffer2& vb, Int row,
462 : Bool forward)
463 : {
464 0 : if(changed(vb, row)) update(vb, row);
465 0 : hasBeenApplied=true;
466 : // now lastUpdateIndex?_p are valid
467 :
468 0 : if (lastUpdateIndex1_p!=lastUpdateIndex2_p)
469 0 : throw(AipsError("BeamSkyJones::apply(SkyComp...) - can only treat homogeneous PB case"));
470 : else {
471 : // for debug
472 : // cout<<endl<<"BeamSkyJones::apply(SkyComp...) index="<<lastUpdateIndex1_p<<" feed="<<vb.feed1()(0)<<" direction=";
473 : // printDirection(cout,lastDirections_p[lastUpdateIndex1_p]); cout<<endl<<endl;
474 : //
475 0 : PBMath myPBMath;
476 0 : MDirection compdir=in.shape().refDirection();
477 0 : MDirection::Types dirType=MDirection::castType(compdir.getRef().getType());
478 0 : if (getPBMath(lastUpdateIndex1_p, myPBMath))
479 0 : return myPBMath.applyPB(in, out,convertDir(vb, lastDirections_p[lastUpdateIndex1_p], dirType) ,
480 0 : Quantity(vb.getFrequency(0,0), "Hz"),
481 0 : lastParallacticAngles_p[lastUpdateIndex1_p],
482 0 : doSquint_p, false, threshold(), forward);
483 : else
484 0 : throw(AipsError("BeamSkyJones::apply(SkyComponent,...) - PBMath not found"));
485 : }
486 : };
487 :
488 :
489 : SkyComponent&
490 0 : BeamSkyJones::applySquare(SkyComponent& in,
491 : SkyComponent& out,
492 : const vi::VisBuffer2& vb, Int row)
493 : {
494 0 : if(changed(vb, row)) update(vb, row);
495 0 : hasBeenApplied=true;
496 : // now lastUpdateIndex?_p are valid
497 :
498 0 : if (lastUpdateIndex1_p!=lastUpdateIndex2_p)
499 0 : throw(AipsError("BeamSkyJones::applySquare(SkyComponent,...) - can only treat homogeneous PB case"));
500 : else {
501 0 : PBMath myPBMath;
502 0 : MDirection compdir=in.shape().refDirection();
503 0 : MDirection::Types dirType=MDirection::castType(compdir.getRef().getType());
504 0 : if (getPBMath(lastUpdateIndex1_p, myPBMath))
505 0 : return myPBMath.applyPB2(in, out, convertDir(vb, lastDirections_p[lastUpdateIndex1_p], dirType),
506 0 : Quantity(vb.getFrequency(0,0), "Hz"),
507 0 : lastParallacticAngles_p[lastUpdateIndex1_p],
508 0 : doSquint_p);
509 : else
510 0 : throw(AipsError("BeamSkyJones::applySquare(SkyComponent,...) - PBMath not found"));
511 : }
512 : };
513 :
514 : // Apply gradient
515 : ImageInterface<Complex>&
516 0 : BeamSkyJones::applyGradient(ImageInterface<Complex>& result,
517 : const vi::VisBuffer2&,
518 : Int)
519 : {
520 0 : return result;
521 : };
522 :
523 :
524 0 : void BeamSkyJones::initializeGradients()
525 : {
526 0 : };
527 :
528 0 : void BeamSkyJones::finalizeGradients()
529 : {
530 0 : };
531 :
532 :
533 : SkyComponent&
534 0 : BeamSkyJones::applyGradient(SkyComponent& result, const vi::VisBuffer2&,
535 : Int)
536 : {
537 0 : return result;
538 : };
539 :
540 0 : void BeamSkyJones::addGradients(const vi::VisBuffer2&, Int,
541 : const Float,
542 : const Float,
543 : const Matrix<Complex>&,
544 : const Matrix<Float>&)
545 0 : {};
546 :
547 : // Solve
548 : /*Bool BeamSkyJones::solve (SkyEquation& se)
549 : {
550 : // Keep compiler quiet
551 : if(&se) {};
552 : return false;
553 : };
554 : */
555 : // return index of compareTelescope, compareAntenna and compareFeed in
556 : // myTelescopes_p; -1 if not found
557 : // if compareAntenna or compareTelescope is -1, this means that the
558 : // PBMath class is the same for all antennae/feeds. An exception will
559 : // be raised, if separate PBMath objects have been assigned by setPBMath
560 : // for different feeds/antennae but -1 is used for query.
561 : //
562 : // It would be good to rename this function to indexBeams as this name
563 : // is more appropriate.
564 : //
565 605 : Int BeamSkyJones::indexTelescope(const String &compareTelescope,
566 : const Int &compareAntenna, const Int &compareFeed) const
567 : {
568 : // for debugging
569 : //cout<<endl<<"BeamSkyJones::indexTelescope tel="<<compareTelescope<<" ant="<<compareAntenna<<" feed="<<compareFeed<<endl<<endl;
570 : //cout<<"Currently "<<myTelescopes_p.nelements()<<" models have been set"<<endl;
571 : //for (uInt i=0; i<myTelescopes_p.nelements(); ++i)
572 : // cout<<i<<" telescope: "<<myTelescopes_p[i]<<" ant:" <<
573 : // myAntennaIDs_p[i]<<" feed: "<<myFeedIDs_p[i]<<endl;
574 : //
575 :
576 605 : DebugAssert(myTelescopes_p.nelements()==myAntennaIDs_p.nelements(),
577 : AipsError);
578 605 : DebugAssert(myTelescopes_p.nelements()==myFeedIDs_p.nelements(),
579 : AipsError);
580 605 : for (uInt i=0; i<myTelescopes_p.nelements(); ++i)
581 192 : if (myTelescopes_p[i] == compareTelescope) {
582 192 : if (compareAntenna==myAntennaIDs_p[i] &&
583 0 : compareFeed==myFeedIDs_p[i]) return i; // -1 will also work
584 192 : if (compareAntenna==myAntennaIDs_p[i]) {
585 0 : if (compareFeed==-1)
586 : throw AipsError("BeamSkyJones::indexTelescope: separate beam models"
587 0 : "have been set up for different feeds and a common model is requested");
588 0 : if (myFeedIDs_p[i]==-1) return i; // valid for all feeds and a given antenna
589 : }
590 192 : if (compareFeed==myFeedIDs_p[i]) {
591 0 : if (compareAntenna==-1)
592 : throw AipsError("BeamSkyJones::indexTelescope: separate beam models"
593 0 : "have been set up for different antennae and a common model is requested");
594 0 : if (myAntennaIDs_p[i]==-1) return i; // valid for all antennae and a given feed
595 : }
596 192 : if (myFeedIDs_p[i]==-1 && myAntennaIDs_p[i]==-1)
597 192 : return i; // valid for all antennae and feeds
598 :
599 : }
600 413 : return -1;
601 : };
602 :
603 : // get the PBMath object; returns false if that one doesn't exist,
604 : // true if it does exist and is OK
605 : // antennaID and feedID default to -1 to preserve the old interface
606 : // TODO: change the interface and make antennaID and feedID the
607 : // second and third parameter respectively to have a better looking code
608 :
609 0 : Bool BeamSkyJones::getPBMath(const String &telescope, PBMath &myPBMath,
610 : const Int &antennaID, const Int &feedID) const
611 : {
612 0 : Int indTel = indexTelescope(telescope,antennaID,feedID);
613 0 : if (indTel >= 0)
614 0 : return getPBMath((uInt)indTel, myPBMath);
615 : else
616 0 : return false; // PBMath not found for this telescope/antenna/feed combination
617 :
618 : };
619 :
620 96 : Bool BeamSkyJones::getPBMath(uInt whichAnt, PBMath &myPBMath) const
621 : {
622 96 : if (whichAnt < myPBMaths_p.nelements() && Int(whichAnt)>=0) {
623 96 : if (myPBMaths_p[whichAnt].ok()) {
624 96 : myPBMath = myPBMaths_p[whichAnt];
625 96 : return true;
626 : } else
627 0 : return false; // whichAnt's PBMath found but not valid
628 : } else
629 0 : return false; // whichAnt's PBMath not found
630 :
631 : };
632 :
633 : // set the PB based on telescope name, antennaID and feedID
634 : // If antennaID or feedID is -1, the PBMath object is set for
635 : // all antennae or feeds, respectively. These are the default
636 : // values to retain the previous interface.
637 : //
638 : // Note. It would be nice to change the interface and make antennaID
639 : // and feedID the second and the third parameter, respectively.
640 215 : void BeamSkyJones::setPBMath(const String &telescope, PBMath &myPBMath,
641 : const Int &antennaID, const Int &feedID)
642 : {
643 : // for debug
644 : // cout<<endl<<"BeamSkyJones::setPBMath tel="<<telescope<<" ant="<<antennaID<<" feed="<<feedID<<endl<<endl;
645 : //
646 :
647 215 : DebugAssert(myTelescopes_p.nelements()==myAntennaIDs_p.nelements(),
648 : AipsError);
649 215 : DebugAssert(myTelescopes_p.nelements()==myFeedIDs_p.nelements(),
650 : AipsError);
651 215 : DebugAssert(myTelescopes_p.nelements()==myPBMaths_p.nelements(),
652 : AipsError);
653 :
654 215 : Bool doRemove=false;
655 215 : if (antennaID==-1 || feedID==-1)
656 : // we have to remove PBMaths for individual antennae/feeds, if they
657 : // were assigned earlier
658 215 : for (uInt i=0; i<myTelescopes_p.nelements(); ++i) {
659 0 : if (doRemove) {
660 : // we have to step back because the previous element
661 : // has been removed
662 0 : --i;
663 0 : doRemove=false;
664 0 : DebugAssert(i<myTelescopes_p.nelements(), AipsError);
665 : }
666 0 : if (myTelescopes_p[i] == telescope) {
667 0 : if (myAntennaIDs_p[i]==-1 && antennaID==-1 &&
668 0 : myFeedIDs_p[i]==-1 && feedID==-1)
669 0 : continue; // to speed things up
670 0 : if ((myAntennaIDs_p[i]!=-1) && (antennaID==-1))
671 : {
672 0 : if (myFeedIDs_p[i]!=-1) myAntennaIDs_p[i]=-1;
673 : // now it's valid for all antennae and a given feed
674 : // and will be replaced later
675 0 : else doRemove=true;
676 : }
677 0 : if ((myFeedIDs_p[i]!=-1) && (feedID==-1))
678 : {
679 0 : if (myAntennaIDs_p[i]!=-1) myFeedIDs_p[i]=-1;
680 : // now it's valid for all feeds at a given antenna
681 : // and will be replaced later
682 0 : else doRemove=true;
683 : }
684 0 : if (doRemove) {
685 0 : myTelescopes_p.remove(i,false);
686 0 : myAntennaIDs_p.remove(i,false);
687 0 : myFeedIDs_p.remove(i,false);
688 0 : myPBMaths_p.remove(i,false);
689 0 : if (lastParallacticAngles_p.nelements())
690 0 : lastParallacticAngles_p.remove(i,false);
691 0 : if (lastDirections_p.nelements())
692 0 : lastDirections_p.remove(i,false);
693 : }
694 : }
695 : }
696 : // there should be no exception at this stage because all bad elements
697 : // should be removed by the code above
698 215 : Int ind = indexTelescope(telescope,antennaID,feedID);
699 215 : if (ind < 0) {
700 215 : ind = myPBMaths_p.nelements();
701 215 : myPBMaths_p.resize(ind+1);
702 215 : myTelescopes_p.resize(ind+1);
703 215 : myTelescopes_p[ind] = telescope;
704 215 : myAntennaIDs_p.resize(ind+1);
705 215 : myAntennaIDs_p[ind] = antennaID;
706 215 : myFeedIDs_p.resize(ind+1);
707 215 : myFeedIDs_p[ind] = feedID;
708 215 : if (lastParallacticAngles_p.nelements())
709 0 : lastParallacticAngles_p.resize(ind+1);
710 215 : if (lastDirections_p.nelements())
711 0 : lastDirections_p.resize(ind+1);
712 : }
713 215 : myPBMaths_p[ind] = myPBMath;
714 215 : if (lastParallacticAngles_p.nelements())
715 0 : lastParallacticAngles_p[ind]=1000.; // to force
716 : // recalculation (it is >> 2pi)
717 215 : };
718 78 : MDirection BeamSkyJones::convertDir(const vi::VisBuffer2& vb, const MDirection& inDir, const MDirection::Types outType){
719 :
720 :
721 78 : if(MDirection::castType(inDir.getRef().getType())==outType){
722 0 : return inDir;
723 : }
724 156 : MPosition pos;
725 156 : String tel("");
726 156 : MSColumns msc(vb.ms());
727 78 : if (vb.subtableColumns().observation().nrow() > 0) {
728 78 : tel = vb.subtableColumns().observation().telescopeName()(vb.observationId()(0));
729 : }
730 138 : if (tel.length() == 0 || !tel.contains("VLA") ||
731 60 : !MeasTable::Observatory(pos,tel)) {
732 : // unknown observatory, use first antenna
733 18 : Int ant1=vb.antenna1()(0);
734 18 : pos=vb.subtableColumns().antenna().positionMeas()(ant1);
735 : }
736 78 : MEpoch::Types timeMType=MEpoch::castType(msc.timeMeas()(0).getRef().getType());
737 156 : Unit timeUnit=Unit(msc.timeMeas().measDesc().getUnits()(0).getName());
738 156 : MEpoch timenow(Quantity(msc.time()(0), timeUnit), timeMType);
739 156 : MeasFrame mFrame(timenow, pos);
740 78 : MDirection::Ref elRef(outType, mFrame);
741 78 : return MDirection::Convert(inDir, elRef)();
742 : }
743 :
744 :
745 :
746 0 : Bool BeamSkyJones::isHomogeneous() const
747 : {
748 : // Hogwash! our "myPBMath_p/myTelescope_p scheme only deals
749 : // with homogeneous pointings. Need to fix this!
750 : // Wait for MS-II
751 0 : return true;
752 : };
753 :
754 :
755 : ImageRegion*
756 0 : BeamSkyJones::extent (const ImageInterface<Complex>& im,
757 : const vi::VisBuffer2& vb,
758 : const Int row,
759 : const Float fPad,
760 : const Int iChan,
761 : const casa::SkyJones::SizeType sizeType)
762 : {
763 0 : if(changed(vb, row)) update(vb, row);
764 0 : DebugAssert(lastUpdateIndex1_p>=0,AipsError); // should be checked in changed/update
765 :
766 0 : PBMath myPBMath;
767 0 : if (getPBMath(lastUpdateIndex1_p, myPBMath)) {
768 0 : return myPBMath.extent(im, lastDirections_p[lastUpdateIndex1_p], row, fPad, iChan, sizeType);
769 : } else {
770 0 : throw(AipsError("BeamSkyJones::extent() - PBMath not found"));
771 : }
772 : };
773 :
774 : ImageRegion*
775 0 : BeamSkyJones::extent (const ImageInterface<Float>& im,
776 : const vi::VisBuffer2& vb,
777 : const Int row,
778 : const Float fPad,
779 : const Int iChan,
780 : const casa::SkyJones::SizeType sizeType)
781 : {
782 0 : if(changed(vb, row)) update(vb, row);
783 0 : DebugAssert(lastUpdateIndex1_p>=0,AipsError); // should be checked in changed/update
784 :
785 0 : PBMath myPBMath;
786 0 : if (getPBMath(lastUpdateIndex1_p, myPBMath)) {
787 0 : return myPBMath.extent(im, lastDirections_p[lastUpdateIndex1_p], row, fPad, iChan, sizeType);
788 : } else {
789 0 : throw(AipsError("BeamSkyJones::extent() - PBMath not found"));
790 : }
791 : };
792 :
793 18 : Int BeamSkyJones::support(const vi::VisBuffer2& vb, const casacore::CoordinateSystem& cs){
794 36 : PBMath myPBMath;
795 :
796 18 : if(changed(vb, 0)) update(vb, 0);
797 18 : if (getPBMath(lastUpdateIndex1_p, myPBMath)) {
798 36 : return myPBMath.support(cs);
799 : } else {
800 0 : throw(AipsError("BeamSkyJones::support() - PBMath not found"));
801 : }
802 :
803 : }
804 :
805 0 : void BeamSkyJones::summary(Int n)
806 : {
807 0 : uInt nPBs = myPBMaths_p.nelements();
808 0 : LogIO os(LogOrigin("BeamSkyJones", "summary"));
809 0 : os << "Beam Summary: "<< LogIO::POST;
810 0 : for (uInt i=0; i< nPBs; ++i) {
811 0 : String pbclass;
812 0 : myPBMaths_p[i].namePBClass(pbclass);
813 0 : os << "Model " << i+1 << " for " << myTelescopes_p[i] <<" ant="<<
814 0 : myAntennaIDs_p[i]<<" feed="<<myFeedIDs_p[i]<< " uses PB: "
815 0 : << pbclass << LogIO::POST;
816 0 : if (n >=0) {
817 0 : myPBMaths_p[i].summary(n);
818 : }
819 : }
820 0 : };
821 :
822 : } //# end of namespace refim
823 :
824 : } //# end of namespace casa
|
