Line data Source code
1 : // -*- C++ -*-
2 : //# PointingOffsets.cc: Implementation of the PointingOffsets class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : //
29 :
30 : #include <msvis/MSVis/VisibilityIterator2.h>
31 : #include <casacore/measures/Measures/MeasTable.h>
32 : #include <casacore/ms/MeasurementSets/MSColumns.h>
33 : #include <synthesis/TransformMachines2/PointingOffsets.h>
34 : #include <synthesis/TransformMachines/SynthesisError.h>
35 : // #include <casa/Logging/LogIO.h>
36 : // #include <casa/Logging/LogSink.h>
37 : // #include <casa/Logging/LogOrigin.h>
38 :
39 : using namespace casacore;
40 : namespace casa{
41 : using namespace vi;
42 : using namespace refim;
43 : //
44 : //----------------------------------------------------------------------
45 : //
46 0 : PointingOffsets& PointingOffsets::operator=(const PointingOffsets& other)
47 : {
48 0 : if(this!=&other)
49 : {
50 0 : imageDC_p = other.imageDC_p;
51 0 : imageObsInfo_p = other.imageObsInfo_p;
52 0 : cachedPointingOffsets_p = other.cachedPointingOffsets_p;
53 0 : cachedAntGridPointingOffsets_p = other.cachedAntGridPointingOffsets_p;
54 : }
55 0 : return *this;
56 : }
57 : //
58 : //----------------------------------------------------------------------
59 : //
60 0 : Vector<Vector<double> >PointingOffsets::findMosaicPointingOffset(const ImageInterface<Complex>& image,
61 : const VisBuffer2& vb, const Bool& doPointing)
62 : {
63 0 : storeImageParams(image,vb);
64 : //where in the image in pixels is this pointing
65 0 : pixFieldGrad_p.resize(2);
66 : int antId,numRow;
67 0 : antId = 0;
68 0 : numRow = 0;
69 0 : Vector<Vector<double> >pixFieldGrad_l;
70 0 : pixFieldGrad_l.resize(1);
71 0 : MDirection dir = vbUtils_p.getPointingDir(vb,antId,numRow,dc_p.directionType(),doPointing);
72 0 : thePix_p = toPix(vb, dir, dir);
73 :
74 0 : pixFieldGrad_p = gradPerPixel(thePix_p);
75 0 : pixFieldGrad_l(0)=(pixFieldGrad_p);
76 :
77 0 : return pixFieldGrad_l;
78 : };
79 : //
80 : //----------------------------------------------------------------------
81 : //
82 0 : Vector<Vector<double> > PointingOffsets::findAntennaPointingOffset(const ImageInterface<Complex>& image,
83 : const vi::VisBuffer2& vb, const Bool& doPointing)
84 : {
85 0 : Vector<Vector<double> >antOffsets;
86 0 : storeImageParams(image,vb);
87 :
88 0 : if (epJ_p.null())
89 : {
90 0 : int numRow_p = vb.nRows();
91 0 : antOffsets.resize(numRow_p); // The array is resized to fit for a given vb
92 0 : if (PO_DEBUG_P==1)
93 0 : cerr << "-------------------------------------------------------------------------" << endl;
94 0 : for (int irow=0; irow<numRow_p;irow++)
95 : {
96 0 : MDirection antDir1 =vbUtils_p.getPointingDir(vb, vb.antenna1()[irow], irow, dc_p.directionType(), doPointing);
97 0 : MDirection antDir2 =vbUtils_p.getPointingDir(vb, vb.antenna2()[irow], irow, dc_p.directionType(), doPointing);
98 :
99 : // MVDirection vbdir=vb.direction1()(0).getValue();
100 0 : casacore::Vector<double> thePixDir1_l, thePixDir2_l;
101 0 : thePixDir1_l = toPix(vb, antDir1, antDir2);
102 0 : thePixDir2_l = toPix(vb, antDir2, antDir1);
103 0 : thePix_p = (thePixDir1_l + thePixDir2_l)/2.;
104 0 : antOffsets(irow) = gradPerPixel(thePix_p);
105 0 : if (PO_DEBUG_P==1)
106 : {
107 0 : cerr << irow << " "
108 0 : << vb.antenna1()[irow] << " " << vb.antenna2()[irow] << " "
109 0 : << antDir1 << " " << antDir2 << " "
110 0 : << vb.direction1()(irow) << " " << vb.direction2()(irow) << " "
111 0 : << toPix(vb, antDir1, antDir2) << " " << toPix(vb, antDir2,antDir1)
112 0 : << endl;
113 :
114 : }
115 : //cerr << irow << " " << antOffsets[irow]<< endl;
116 : }
117 0 : if (PO_DEBUG_P==1)
118 0 : cerr << "=========================================================================" << endl;
119 : }
120 : else
121 : {
122 0 : throw(SynthesisFTMachineError("PointingOffsets::findAntennaPointingOffset(): Antenna pointing CalTable support not yet implemented"));
123 : }
124 0 : return antOffsets;
125 : }
126 : //
127 : //----------------------------------------------------------------------
128 : //
129 0 : void PointingOffsets::fetchPointingOffset(const ImageInterface<Complex>& image,
130 : const VisBuffer2& vb, const Bool doPointing)
131 : {
132 0 : setDoPointing(doPointing);
133 0 : if (!doPointing)
134 : {
135 0 : cachedPointingOffsets_p.assign(findMosaicPointingOffset(image,vb,doPointing));
136 : }
137 : else
138 : {
139 0 : cachedPointingOffsets_p.assign(findAntennaPointingOffset(image,vb,doPointing));
140 : }
141 0 : }
142 : //
143 : //----------------------------------------------------------------------
144 : //
145 0 : Vector<double> PointingOffsets::gradPerPixel(const Vector<double>& p)
146 : {
147 0 : Vector<double> gPP(2);
148 0 : gPP(0) = p(0) - double(nx_p/2);
149 0 : gPP(1) = p(1) - double(ny_p/2);
150 :
151 0 : gPP(0) = -gPP(0)*2.0*C::pi/double(nx_p)/double(convOversampling_p);
152 0 : gPP(1) = -gPP(1)*2.0*C::pi/double(ny_p)/double(convOversampling_p);
153 :
154 0 : return gPP;
155 : }
156 : //
157 : //----------------------------------------------------------------------
158 : //
159 0 : Vector<double>& PointingOffsets::toPix(const VisBuffer2& vb,
160 : const MDirection& dir1,
161 : const MDirection& dir2)
162 : {
163 0 : thePix_p.resize(2);
164 0 : if(dc_p.directionType() != MDirection::castType(dir1.getRef().getType()))
165 : {
166 : //pointToPix_p.setModel(theDir);
167 :
168 0 : MEpoch timenow(Quantity(vb.time()(0), timeUnit_p), timeMType_p);
169 0 : pointFrame_p.resetEpoch(timenow);
170 : //////////////////////////
171 : //pointToPix holds pointFrame_p by reference...
172 : //thus good to go for conversion
173 0 : direction1_p=pointToPix_p(dir1);
174 0 : direction2_p=pointToPix_p(dir2);
175 0 : dc_p.toPixel(thePix_p, direction1_p);
176 : //cerr<<" thePix_P from one: "<<thePix_p<< " " << dir1.getRef().getType()<<endl;
177 :
178 : }
179 : else
180 : {
181 0 : direction1_p=dir1;
182 0 : direction2_p=dir2;
183 0 : dc_p.toPixel(thePix_p, dir1);
184 : //cerr<<" thePix_P from two: "<<thePix_p<< " " << dir1.getRef().getType()<<endl;
185 : }
186 : // Return the internal variable, just to make code more readable
187 : // at the place of the call.
188 0 : return thePix_p;
189 : };
190 :
191 : //
192 : //----------------------------------------------------------------------
193 : //
194 :
195 0 : vector<vector<double> > PointingOffsets::fetchAntOffsetToPix(const VisBuffer2& vb, const Bool doPointing)
196 : {
197 : // Int numRow_p = vb.nRows();
198 0 : vector<vector<double> > pix_l;
199 :
200 0 : vector<int> ant1, ant2;
201 0 : ant1 = (vb.antenna1()).tovector();
202 0 : ant2 = (vb.antenna2()).tovector();
203 0 : ant1.insert(ant1.end(),ant2.begin(),ant2.end());
204 0 : sort(ant1.begin(),ant1.end());
205 0 : auto itr = unique(ant1.begin(),ant1.end());
206 0 : ant1.resize(distance(ant1.begin(),itr));
207 :
208 : // cerr <<"ant1.size()" << ant1.size() << endl;
209 0 : pix_l.resize(2);
210 0 : pix_l[0].resize(ant1.size(),0);
211 0 : pix_l[1].resize(ant1.size(),0);
212 :
213 0 : MVDirection vbdir=vb.direction1()(0).getValue();
214 0 : for (unsigned int nant=0; nant< ant1.size();nant++)
215 : {
216 :
217 0 : MDirection antDir1 =vbUtils_p.getPointingDir(vb, ant1[nant], 0, dc_p.directionType(), doPointing);
218 : // MDirection antDir2 =vbUtils_p.getPointingDir(vb, vb.antenna2()[irow], irow, dc_p.directionType(), doPointing);
219 0 : Vector<double> tmp = toPix(vb, antDir1, vbdir);
220 0 : pix_l[0][nant]=tmp[0];
221 0 : pix_l[1][nant]=tmp[1];
222 : // cerr<< "Ant : "<< ant1[nant]<< " Offsets : "<< pix_l[0][nant] << " " << pix_l[1][nant]<<endl;
223 : // tmp = toPix(vb, antDir2, vbdir);
224 : // pix_l[2][irow]=tmp[0];
225 : // pix_l[3][irow]=tmp[1];
226 : }
227 0 : return pix_l;
228 : }
229 :
230 :
231 : //
232 : //----------------------------------------------------------------------
233 : //
234 0 : void PointingOffsets::storeImageParams(const ImageInterface<Complex>& iimage,
235 : const VisBuffer2& vb)
236 : {
237 : //image signature changed...rather simplistic for now
238 0 : if((iimage.shape().product() != nx_p*ny_p*nchan_p*npol_p) || nchan_p < 1)
239 : {
240 0 : csys_p=iimage.coordinates();
241 0 : int coordIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
242 0 : AlwaysAssert(coordIndex>=0, AipsError);
243 0 : directionIndex_p=coordIndex;
244 0 : dc_p=csys_p.directionCoordinate(directionIndex_p);
245 0 : ObsInfo imInfo=csys_p.obsInfo();
246 0 : String tel= imInfo.telescope();
247 0 : MPosition pos;
248 0 : ROMSColumns mscol(vb.ms());
249 0 : if (vb.subtableColumns().observation().nrow() > 0)
250 : {
251 0 : tel = vb.subtableColumns().observation().telescopeName()
252 0 : (mscol.observationId()(0));
253 : }
254 0 : if (tel.length() == 0 || !tel.contains("VLA") ||
255 0 : !MeasTable::Observatory(pos,tel))
256 : {
257 : // unknown observatory, use first antenna
258 0 : int ant1 = vb.antenna1()(0);
259 0 : pos=vb.subtableColumns().antenna().positionMeas()(ant1);
260 : }
261 : //cout << "TELESCOPE " << tel << endl;
262 : //Store this to build epochs via the time access of visbuffer later
263 :
264 0 : timeMType_p=MEpoch::castType(mscol.timeMeas()(0).getRef().getType());
265 0 : timeUnit_p=Unit(mscol.timeMeas().measDesc().getUnits()(0).getName());
266 : // timeUnit_p=Unit("s");
267 : //cout << "UNIT " << timeUnit_p.getValue() << " name " << timeUnit_p.getName() << endl;
268 0 : pointFrame_p=MeasFrame(imInfo.obsDate(), pos);
269 0 : MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
270 : //For now we set the conversion from this direction
271 0 : pointToPix_p=MDirection::Convert( MDirection(), elRef);
272 0 : nx_p=iimage.shape()(coordIndex);
273 0 : ny_p=iimage.shape()(coordIndex+1);
274 0 : coordIndex=csys_p.findCoordinate(Coordinate::SPECTRAL);
275 0 : int pixAxis=csys_p.pixelAxes(coordIndex)[0];
276 0 : nchan_p=iimage.shape()(pixAxis);
277 0 : coordIndex=csys_p.findCoordinate(Coordinate::STOKES);
278 0 : pixAxis=csys_p.pixelAxes(coordIndex)[0];
279 0 : npol_p=iimage.shape()(pixAxis);
280 : }
281 0 : }
282 : };
|
