Line data Source code
1 : // -*- C++ -*-
2 : //# EVLAAperture.cc: Implementation of the EVLAAperture 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 : #include <msvis/MSVis/VisBuffer2.h>
30 : #include <msvis/MSVis/VisibilityIterator2.h>
31 : #include <synthesis/TransformMachines2/EVLAAperture.h>
32 : #include <synthesis/TransformMachines/SynthesisError.h>
33 : #include <synthesis/TransformMachines/BeamCalc.h>
34 : #include <synthesis/TransformMachines2/WTerm.h>
35 : #include <synthesis/TransformMachines2/VLACalcIlluminationConvFunc.h>
36 : #include <synthesis/TransformMachines2/Utils.h>
37 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/ms/MeasurementSets/MSColumns.h>
41 : //
42 : //---------------------------------------------------------------------
43 : //
44 :
45 : using namespace casacore;
46 : namespace casa{
47 : using namespace vi;
48 : using namespace refim;
49 : using namespace SynthesisUtils;
50 0 : EVLAAperture& EVLAAperture::operator=(const EVLAAperture& other)
51 : {
52 0 : if(this!=&other)
53 : {
54 : // ConvolutionFunction::operator=(other);
55 0 : logIO_p = other.logIO_p;
56 : // setParams(other.polMap_p_base, other.feedStokes_p);
57 0 : setPolMap(other.polMap_p_base);
58 0 : Diameter_p=other.Diameter_p;
59 0 : Nant_p=other.Nant_p;
60 0 : HPBW=other.HPBW;
61 0 : sigma=other.sigma;
62 : }
63 0 : return *this;
64 : }
65 0 : Int EVLAAperture::getVLABandID(Double& vbRefFreq,String&telescopeName, const CoordinateSystem& skyCoord)
66 : {
67 0 : LogIO log_l(LogOrigin("EVLAAperture", "getVLABandID[R&D]"));
68 :
69 0 : Double refFreq = skyCoord.spectralCoordinate(skyCoord.findCoordinate(Coordinate::SPECTRAL)).referenceValue()(0);
70 0 : cerr << "getVLABand (Global VB Ref. min, CF Ref.): " << vbRefFreq << " , " << refFreq << endl;
71 0 : if (telescopeName=="VLA")
72 : {
73 0 : if ((refFreq >=1.34E9) && (refFreq <=1.73E9)) return BeamCalc_VLA_L;
74 0 : else if ((refFreq >=4.5E9) && (refFreq <=5.0E9)) return BeamCalc_VLA_C;
75 0 : else if ((refFreq >=8.0E9) && (refFreq <=8.8E9)) return BeamCalc_VLA_X;
76 0 : else if ((refFreq >=14.4E9) && (refFreq <=15.4E9)) return BeamCalc_VLA_U;
77 0 : else if ((refFreq >=22.0E9) && (refFreq <=24.0E9)) return BeamCalc_VLA_K;
78 0 : else if ((refFreq >=40.0E9) && (refFreq <=50.0E9)) return BeamCalc_VLA_Q;
79 0 : else if ((refFreq >=100E6) && (refFreq <=300E6)) return BeamCalc_VLA_4;
80 : }
81 : else
82 0 : if (telescopeName=="EVLA")
83 : {
84 0 : if ((refFreq >= 0.9E9) && (refFreq <= 2.1E9)) return BeamCalc_EVLA_L;
85 0 : else if ((refFreq >=2.0E9) && (refFreq <=4.0E9)) return BeamCalc_EVLA_S;
86 0 : else if ((refFreq >=4.0E9) && (refFreq <=8.0E9)) return BeamCalc_EVLA_C;
87 0 : else if ((refFreq >=8.0E9) && (refFreq <=12.0E9)) return BeamCalc_EVLA_X;
88 0 : else if ((refFreq >=12.0E9) && (refFreq <=18.0E9)) return BeamCalc_EVLA_U;
89 0 : else if ((refFreq >=18.0E9) && (refFreq <=26.5E9)) return BeamCalc_EVLA_K;
90 0 : else if ((refFreq >=26.5E9) && (refFreq <=40.8E9)) return BeamCalc_EVLA_A;
91 0 : else if ((refFreq >=40.0E9) && (refFreq <=50.0E9)) return BeamCalc_EVLA_Q;
92 : }
93 0 : ostringstream mesg;
94 0 : log_l << telescopeName << "/" << refFreq << "(Hz) combination not recognized." << LogIO::EXCEPTION;
95 0 : return -1;
96 : }
97 :
98 0 : void EVLAAperture::cacheVBInfo(const String& telescopeName,
99 : const Float& diameter)
100 : {
101 0 : telescopeName_p=telescopeName;
102 0 : Diameter_p=diameter;
103 0 : }
104 :
105 0 : void EVLAAperture::cacheVBInfo(const VisBuffer2& vb)
106 : {
107 0 : const Vector<String> telescopeNames=vb.subtableColumns()
108 0 : .observation().telescopeName().getColumn();
109 0 : for(uInt nt=0;nt<telescopeNames.nelements();nt++)
110 : {
111 0 : if ((telescopeNames(nt) != "VLA") && (telescopeNames(nt) != "EVLA"))
112 : {
113 0 : String mesg="We can handle only (E)VLA antennas for now.\n";
114 0 : mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
115 0 : SynthesisError err(mesg);
116 0 : throw(err);
117 : }
118 0 : if (telescopeNames(nt) != telescopeNames(0))
119 : {
120 0 : String mesg="We do not (yet) handle multiple telescopes for A-Projection!\n";
121 0 : mesg += "Not yet a \"priority\"!!";
122 0 : SynthesisError err(mesg);
123 0 : throw(err);
124 : }
125 : }
126 0 : telescopeName_p=telescopeNames[0];
127 :
128 : // MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
129 : //Freq = vb.msColumns().spectralWindow().refFrequency()(0);
130 0 : Diameter_p=0;
131 0 : Nant_p = vb.subtableColumns().antenna().nrow();
132 0 : for (Int i=0; i < Nant_p; i++)
133 0 : if (!vb.subtableColumns().antenna().flagRow()(i))
134 : {
135 0 : Diameter_p = vb.subtableColumns().antenna().dishDiameter().getColumn()(i);
136 0 : break;
137 : }
138 0 : if (Diameter_p == 0)
139 : {
140 0 : logIO() << LogOrigin("EVLAAperture", "cacheVBInfo")
141 : << "No valid or finite sized antenna found in the antenna table. "
142 : << "Assuming diameter = 25m."
143 : << LogIO::WARN
144 0 : << LogIO::POST;
145 0 : Diameter_p=25.0;
146 : }
147 0 : cacheVBInfo(telescopeNames[0], Diameter_p);
148 0 : }
149 :
150 0 : Int EVLAAperture::getBandID(const Double& freq, const String& telescopeName, const String& bandName)
151 : {
152 0 : Int bandID=0;
153 0 : if (!isNoOp())
154 : {
155 : // First #-separated token in bandName_p is the name of the band used
156 0 : Vector<String> tokens = SynthesisUtils::parseBandName(bandName);
157 0 : bandID = BeamCalc::Instance()->getBandID(freq,telescopeName,tokens(0));
158 : }
159 :
160 0 : return bandID;
161 : };
162 :
163 0 : int EVLAAperture::getVisParams(const VisBuffer2& vb,const CoordinateSystem& /*im*/)
164 : {
165 0 : throw(AipsError("EVLAAperture::getVisParams() called"));
166 : Double Freq;
167 : cacheVBInfo(vb);
168 :
169 : Freq = vb.subtableColumns().spectralWindow().refFrequency()(0);
170 : Double Lambda=C::c/Freq;
171 : HPBW = Lambda/(Diameter_p*sqrt(log(2.0)));
172 : sigma = 1.0/(HPBW*HPBW);
173 :
174 : return getBandID(Freq, telescopeName_p,bandName_p);
175 : // Int bandID=0;
176 : // if (!isNoOp())
177 : // bandID = BeamCalc::Instance()->getBandID(Freq,telescopeName_p);
178 :
179 : // return bandID;
180 : }
181 :
182 0 : Int EVLAAperture::makePBPolnCoords(const VisBuffer2&vb,
183 : const Int& convSize,
184 : const Int& convSampling,
185 : const CoordinateSystem& skyCoord,
186 : const Int& skyNx, const Int& /*skyNy*/,
187 : CoordinateSystem& feedCoord)
188 : // Vector<Int>& cfStokes)
189 : {
190 0 : feedCoord = skyCoord;
191 : //
192 : // Make a two dimensional image to calculate auto-correlation of
193 : // the ideal illumination pattern. We want this on a fine grid in
194 : // the UV plane
195 : //
196 0 : Int directionIndex=skyCoord.findCoordinate(Coordinate::DIRECTION);
197 : // cout<<"Direction index is "<< directionIndex<<endl;
198 0 : AlwaysAssert(directionIndex>=0, AipsError);
199 0 : DirectionCoordinate dc=skyCoord.directionCoordinate(directionIndex);
200 0 : Vector<Double> sampling;
201 0 : sampling = dc.increment();
202 : // cout<<"Image sampling set to : "<<sampling<<endl;
203 0 : sampling*=Double(convSampling);
204 0 : sampling*=Double(skyNx)/Double(convSize);
205 0 : dc.setIncrement(sampling);
206 : // cout<<"Resized sampling is : "<<sampling<<endl;
207 :
208 0 : Vector<Double> unitVec(2);
209 0 : unitVec=convSize/2;
210 0 : dc.setReferencePixel(unitVec);
211 :
212 : // Set the reference value to that of the image
213 0 : feedCoord.replaceCoordinate(dc, directionIndex);
214 :
215 : //
216 : // Make an image with circular polarization axis.
217 : //
218 0 : Int NPol=0,M,N=0;
219 0 : M=polMap_p_base.nelements();
220 0 : for(Int i=0;i<M;i++) if (polMap_p_base(i) > -1) NPol++;
221 0 : Vector<Int> poln(NPol);
222 :
223 : Int index;
224 0 : Vector<Int> inStokes;
225 0 : index = feedCoord.findCoordinate(Coordinate::STOKES);
226 0 : inStokes = feedCoord.stokesCoordinate(index).stokes();
227 0 : N = 0;
228 : try
229 : {
230 : // cerr << "### " << polMap_p_base << " " << vb.corrType() << endl;
231 0 : for(Int i=0;i<M;i++) if (polMap_p_base(i) > -1) {poln(N) = vb.correlationTypes()(i);N++;}
232 0 : StokesCoordinate polnCoord(poln);
233 0 : Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
234 0 : feedCoord.replaceCoordinate(polnCoord,StokesIndex);
235 : // cfStokes = poln;
236 : }
237 0 : catch(AipsError& x)
238 : {
239 : throw(SynthesisFTMachineError("Likely cause: Discrepancy between the poln. "
240 0 : "axis of the data and the image specifications."));
241 : }
242 :
243 0 : return NPol;
244 : }
245 :
246 0 : Bool EVLAAperture::findSupport(Array<Complex>& func, Float& threshold,Int& origin, Int& R)
247 : {
248 : Double NSteps;
249 0 : Int PixInc=1;
250 0 : Vector<Complex> vals;
251 0 : IPosition ndx(4,origin,0,0,0);
252 0 : Bool found=false;
253 0 : IPosition cfShape=func.shape();
254 0 : Int convSize = cfShape(0);
255 0 : for(R=convSize/4;R>1;R--)
256 : {
257 0 : NSteps = 90*R/PixInc; //Check every PixInc pixel along a
258 : //circle of radious R
259 0 : vals.resize((Int)(NSteps+0.5));
260 0 : vals=0;
261 0 : for(Int th=0;th<NSteps;th++)
262 : {
263 0 : ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
264 0 : ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
265 :
266 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
267 0 : vals(th)=func(ndx);
268 : }
269 0 : if (max(abs(vals)) > threshold)
270 0 : {found=true;break;}
271 : }
272 0 : return found;
273 : }
274 :
275 0 : void EVLAAperture::makeFullJones(ImageInterface<Complex>& pbImage,
276 : const VisBuffer2& vb,
277 : Bool doSquint, Int& bandID, Double freqVal)
278 : {
279 :
280 0 : if (!isNoOp())
281 : {
282 0 : VLACalcIlluminationConvFunc vlaPB;
283 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
284 0 : vlaPB.setMaximumCacheSize(cachesize);
285 0 : bandID = getBandID(freqVal,telescopeName_p,bandName_p);
286 0 : vlaPB.makeFullJones(pbImage,vb, doSquint, bandID, freqVal);
287 : }
288 0 : }
289 :
290 0 : void EVLAAperture::applySky(ImageInterface<Complex>& outImages,
291 : const VisBuffer2& vb,
292 : const Bool doSquint,
293 : const Int& cfKey,
294 : const Int& muellerTerm,
295 : const Double freqVal)
296 : {
297 : // (void)cfKey;
298 : // if (!isNoOp())
299 : // {
300 : // VLACalcIlluminationConvFunc vlaPB;
301 : // Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
302 : // vlaPB.setMaximumCacheSize(cachesize);
303 : // bandID = getBandID(freqVal,telescopeName_p);
304 : // // cout<<"EVLAAperture : muellerTerm"<<muellerTerm <<"\n";
305 : // //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
306 : // Double pa=getPA(vb);
307 : // vlaPB.applyPB(outImages, pa, doSquint,bandID,muellerTerm,freqVal);
308 :
309 : // }
310 0 : Double pa=getPA(vb);
311 0 : applySky(outImages,pa,doSquint,cfKey,muellerTerm, freqVal);
312 0 : }
313 : //
314 : // This one without VB. Should become the default (and the only one!).
315 : //
316 0 : void EVLAAperture::applySky(ImageInterface<Complex>& outImages,
317 : const Double& pa,
318 : const Bool doSquint,
319 : const Int& cfKey,
320 : const Int& muellerTerm,
321 : const Double freqVal)
322 : {
323 : (void)cfKey;
324 0 : if (!isNoOp())
325 : {
326 0 : VLACalcIlluminationConvFunc vlaPB;
327 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
328 0 : vlaPB.setMaximumCacheSize(cachesize);
329 : Int bandID;
330 0 : bandID = getBandID(freqVal,telescopeName_p,bandName_p);
331 : //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
332 0 : Double pa_l=pa; // Due to goofup in making sure complier type checking does not come in the way!
333 0 : vlaPB.applyPB(outImages, pa_l, doSquint,bandID,muellerTerm,freqVal);
334 :
335 : }
336 0 : }
337 :
338 0 : void EVLAAperture::applySky(ImageInterface<Float>& outImages,
339 : const VisBuffer2& vb,
340 : const Bool doSquint,
341 : const Int& cfKey,
342 : const Int& muellerTerm,
343 : const Double freqVal)
344 : {
345 : (void)cfKey;
346 : (void)muellerTerm;
347 0 : if (!isNoOp())
348 : {
349 0 : VLACalcIlluminationConvFunc vlaPB;
350 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
351 0 : vlaPB.setMaximumCacheSize(cachesize);
352 : Int bandID;
353 0 : bandID = getBandID(freqVal,telescopeName_p,bandName_p);
354 :
355 0 : Double pa=getPA(vb);
356 0 : vlaPB.applyPB(outImages, pa, bandID, doSquint,freqVal);
357 : }
358 0 : }
359 :
360 : };
|
