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