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: casa-feedback@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 <casacore/coordinates/Coordinates/DirectionCoordinate.h>
30 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
31 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
32 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
33 : #include <casacore/images/Images/SubImage.h>
34 : #include <casacore/ms/MeasurementSets/MSColumns.h>
35 : #include <msvis/MSVis/VisBuffer2.h>
36 : #include <msvis/MSVis/VisibilityIterator2.h>
37 :
38 : #include <synthesis/TransformMachines/SynthesisError.h>
39 : #include <synthesis/TransformMachines2/EVLAAperture.h>
40 : #include <synthesis/TransformMachines2/Utils.h>
41 : #include <synthesis/TransformMachines2/VLACalcIlluminationConvFunc.h>
42 : #include <synthesis/TransformMachines2/WTerm.h>
43 : #include <synthesis/Utilities/FFT2D.h>
44 : //
45 : //---------------------------------------------------------------------
46 : //
47 :
48 : using namespace casacore;
49 : namespace casa {
50 : using namespace vi;
51 : using namespace refim;
52 : using namespace SynthesisUtils;
53 :
54 0 : EVLAAperture::EVLAAperture() : AzElAperture(), polMap_p(), feedStokes_p() {
55 0 : telescopeName_p = "EVLA";
56 0 : Diameter_p = 25.0;
57 0 : }
58 :
59 0 : EVLAAperture &EVLAAperture::operator=(const EVLAAperture &other) {
60 0 : if (this != &other) {
61 : // ConvolutionFunction::operator=(other);
62 0 : logIO_p = other.logIO_p;
63 : // setParams(other.polMap_p_base, other.feedStokes_p);
64 0 : setPolMap(other.polMap_p_base);
65 0 : telescopeName_p = other.telescopeName_p;
66 0 : Diameter_p = other.Diameter_p;
67 0 : Nant_p = other.Nant_p;
68 0 : HPBW = other.HPBW;
69 0 : sigma = other.sigma;
70 : }
71 0 : return *this;
72 : }
73 0 : Int EVLAAperture::getVLABandID(Double &vbRefFreq, String &telescopeName,
74 : const CoordinateSystem &skyCoord) {
75 0 : LogIO log_l(LogOrigin("EVLAAperture", "getVLABandID[R&D]"));
76 :
77 : Double refFreq =
78 0 : skyCoord.spectralCoordinate(skyCoord.findCoordinate(Coordinate::SPECTRAL))
79 0 : .referenceValue()(0);
80 : //cerr << "getVLABand (Global VB Ref. min, CF Ref.): " << vbRefFreq << " , "
81 : // << refFreq << endl;
82 0 : if (telescopeName == "VLA") {
83 0 : if ((refFreq >= 1.34E9) && (refFreq <= 1.73E9))
84 0 : return BeamCalc_VLA_L;
85 0 : else if ((refFreq >= 4.5E9) && (refFreq <= 5.0E9))
86 0 : return BeamCalc_VLA_C;
87 0 : else if ((refFreq >= 8.0E9) && (refFreq <= 8.8E9))
88 0 : return BeamCalc_VLA_X;
89 0 : else if ((refFreq >= 14.4E9) && (refFreq <= 15.4E9))
90 0 : return BeamCalc_VLA_U;
91 0 : else if ((refFreq >= 22.0E9) && (refFreq <= 24.0E9))
92 0 : return BeamCalc_VLA_K;
93 0 : else if ((refFreq >= 40.0E9) && (refFreq <= 50.0E9))
94 0 : return BeamCalc_VLA_Q;
95 0 : else if ((refFreq >= 100E6) && (refFreq <= 300E6))
96 0 : return BeamCalc_VLA_4;
97 0 : } else if (telescopeName == "EVLA") {
98 0 : if ((refFreq >= 0.9E9) && (refFreq <= 2.1E9))
99 0 : return BeamCalc_EVLA_L;
100 0 : else if ((refFreq >= 2.0E9) && (refFreq <= 4.0E9))
101 0 : return BeamCalc_EVLA_S;
102 0 : else if ((refFreq >= 4.0E9) && (refFreq <= 8.0E9))
103 0 : return BeamCalc_EVLA_C;
104 0 : else if ((refFreq >= 8.0E9) && (refFreq <= 12.0E9))
105 0 : return BeamCalc_EVLA_X;
106 0 : else if ((refFreq >= 12.0E9) && (refFreq <= 18.0E9))
107 0 : return BeamCalc_EVLA_U;
108 0 : else if ((refFreq >= 18.0E9) && (refFreq <= 26.5E9))
109 0 : return BeamCalc_EVLA_K;
110 0 : else if ((refFreq >= 26.5E9) && (refFreq <= 40.8E9))
111 0 : return BeamCalc_EVLA_A;
112 0 : else if ((refFreq >= 40.0E9) && (refFreq <= 50.0E9))
113 0 : return BeamCalc_EVLA_Q;
114 : }
115 0 : ostringstream mesg;
116 : log_l << telescopeName << "/" << refFreq << "(Hz) combination not recognized."
117 0 : << LogIO::EXCEPTION;
118 0 : return -1;
119 0 : }
120 :
121 0 : void EVLAAperture::setApertureParams(ApertureCalcParams &ap, const Float &Freq,
122 : const Float &pa, const Int &bandID,
123 : const IPosition &skyShape,
124 : const Vector<Double> &uvIncr) {
125 0 : Double Lambda = C::c / Freq;
126 :
127 0 : ap.oversamp = 3;
128 0 : ap.pa = pa;
129 0 : ap.band = bandID;
130 0 : ap.freq = Freq / 1E9;
131 0 : ap.nx = skyShape(0);
132 0 : ap.ny = skyShape(1);
133 0 : ap.dx = abs(uvIncr(0) * Lambda);
134 0 : ap.dy = abs(uvIncr(1) * Lambda);
135 0 : ap.x0 = -(ap.nx / 2) * ap.dx;
136 0 : ap.y0 = -(ap.ny / 2) * ap.dy;
137 : //cerr << "pa= " << ap.pa << " band " << ap.band << " freq " << ap.freq
138 : // << " nx ny " << ap.nx << " " << ap.ny << " dx dy " << ap.dx << " "
139 : // << ap.dy << endl;
140 : //
141 : // If cross-hand pols. are requested, we need to compute both
142 : // the parallel-hand aperture illuminations.
143 : //
144 : // if ((inStokes == Stokes::RL) || (inStokes == Stokes::LR))
145 : {
146 : // IPosition apShape(ap.aperture->shape());
147 : // cerr << "APshape " << apShape << endl;
148 : // apShape(3)=4;
149 : //cerr << "APSHPE=" << skyShape << endl;
150 :
151 : //ap.aperture->resize(skyShape);
152 : }
153 0 : }
154 0 : void EVLAAperture::cacheVBInfo(const String &telescopeName,
155 : const Float &diameter) {
156 0 : telescopeName_p = telescopeName;
157 0 : Diameter_p = diameter;
158 0 : }
159 :
160 0 : void EVLAAperture::cacheVBInfo(const VisBuffer2 &vb) {
161 : const Vector<String> telescopeNames =
162 0 : vb.subtableColumns().observation().telescopeName().getColumn();
163 0 : for (uInt nt = 0; nt < telescopeNames.nelements(); nt++) {
164 0 : if ((telescopeNames(nt) != "VLA") && (telescopeNames(nt) != "EVLA")) {
165 0 : String mesg = "We can handle only (E)VLA antennas for now.\n";
166 0 : mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
167 0 : SynthesisError err(mesg);
168 0 : throw(err);
169 0 : }
170 0 : if (telescopeNames(nt) != telescopeNames(0)) {
171 : String mesg =
172 0 : "We do not (yet) handle multiple telescopes for A-Projection!\n";
173 0 : mesg += "Not yet a \"priority\"!!";
174 0 : SynthesisError err(mesg);
175 0 : throw(err);
176 0 : }
177 : }
178 0 : telescopeName_p = telescopeNames[0];
179 :
180 : // MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
181 : // Freq = vb.msColumns().spectralWindow().refFrequency()(0);
182 0 : Diameter_p = 0;
183 0 : Nant_p = vb.subtableColumns().antenna().nrow();
184 0 : for (Int i = 0; i < Nant_p; i++)
185 0 : if (!vb.subtableColumns().antenna().flagRow()(i)) {
186 0 : Diameter_p = vb.subtableColumns().antenna().dishDiameter().getColumn()(i);
187 0 : break;
188 : }
189 0 : if (Diameter_p == 0) {
190 0 : logIO() << LogOrigin("EVLAAperture", "cacheVBInfo")
191 : << "No valid or finite sized antenna found in the antenna table. "
192 0 : << "Assuming diameter = 25m." << LogIO::WARN << LogIO::POST;
193 0 : Diameter_p = 25.0;
194 : }
195 0 : cacheVBInfo(telescopeNames[0], Diameter_p);
196 0 : }
197 :
198 0 : Int EVLAAperture::getBandID(const Double &freq,
199 : const String & /*telescopeName*/,
200 : const String &bandName) {
201 0 : Int bandID = 0;
202 0 : if (!isNoOp()) {
203 : // First #-separated token in bandName_p is the name of the band used
204 0 : Vector<String> tokens = SynthesisUtils::parseBandName(bandName);
205 0 : Double elfreq = freq;
206 0 : if(telescopeName_p=="VLA" && freq < 1.34e9) //Some unit test data for VLA are below 1.0 GHz
207 0 : elfreq = 1.4e9;
208 0 : bandID = BeamCalc::Instance()->getBandID(elfreq, telescopeName_p, tokens(0));
209 0 : }
210 :
211 0 : return bandID;
212 : };
213 0 : String EVLAAperture::getVLABandName(const Double& freq, const String& telescopeName) {
214 0 : double tol = FLT_EPSILON;
215 0 : String bandName = "EVLA_L";
216 0 : if (telescopeName == "VLA") {
217 : // if ((freq >= 1.34E9) && (freq <= 1.73E9)) some unit test data goes from 1 to 2 GHz for VLA !
218 :
219 0 : if((freq >= (9E8-tol)) && (freq <= (2.0E9+tol)))
220 :
221 0 : bandName = "VLA_L";
222 0 : else if ((freq >= (4.5E9-tol)) && (freq <= (5.0E9+tol)))
223 0 : bandName = "VLA_C";
224 0 : else if ((freq >= (8.0E9-tol)) && (freq <= (8.8E9+tol)))
225 0 : bandName = "VLA_X";
226 0 : else if ((freq >= (14.4E9-tol)) && (freq <= (15.4E9+tol)))
227 0 : bandName = "VLA_U";
228 0 : else if ((freq >= (22.0E9-tol)) && (freq <= (24.0E9+tol)))
229 0 : bandName = "VLA_K";
230 0 : else if ((freq >= (40.0E9-tol)) && (freq <= (50.0E9+tol)))
231 0 : bandName = "VLA_Q";
232 0 : else if ((freq >= (30E6-tol)) && (freq <= (100E6+tol)))
233 0 : bandName = "VLA_4";
234 : else
235 0 : throw(
236 0 : AipsError("Don't know VLA band for frequency=" + String::toString(freq)));
237 0 : } else if (telescopeName == "EVLA") {
238 0 : if (freq > (9e8-tol) && freq <= 2.0e9)
239 0 : bandName = "EVLA_L";
240 0 : else if (freq > 2.0e9 && freq <= 4.0e9)
241 0 : bandName = "EVLA_S";
242 0 : else if (freq > 4.0e9 && freq <= 8.0e9)
243 0 : bandName = "EVLA_C";
244 0 : else if (freq > 8.e9 && freq <= 12.0e9)
245 0 : bandName = "EVLA_X";
246 0 : else if (freq > 12.0e9 && freq <= 18.0e9)
247 0 : bandName = "EVLA_U";
248 0 : else if (freq > 18.0e9 && freq <= 26.0e9)
249 0 : bandName = "EVLA_K";
250 0 : else if (freq > 26.e9 && freq <= 40.0e9)
251 0 : bandName = "EVLA_A";
252 0 : else if (freq > 40.0e9 && freq <= (50.0e9+tol))
253 0 : bandName = "EVLA_Q";
254 : else
255 0 : throw(
256 0 : AipsError("Don't know EVLA band for frequency=" + String::toString(freq)));
257 : } else {
258 0 : throw(AipsError("Don't know telescope " + telescopeName));
259 : }
260 0 : return bandName;
261 0 : }
262 0 : Int EVLAAperture::getBandID(const Double &freq, const String& bandName) {
263 0 : Int bandID = 0;
264 0 : if(bandName=="")
265 0 : bandName_p = getVLABandName(freq, telescopeName_p);
266 : else
267 0 : bandName_p=bandName;
268 :
269 0 : if (!isNoOp()) {
270 : // First #-separated token in bandName_p is the name of the band used
271 : // Vector<String> tokens = SynthesisUtils::parseBandName(bandName_p);
272 : // cerr << "TOKENS " << tokens << endl;
273 0 : Double elfreq = freq;
274 0 : if (telescopeName_p == "VLA" &&
275 0 : freq < 1.34e9) // Some unit test data for VLA are below 1.0 GHz
276 0 : elfreq = 1.4e9;
277 0 : bandID = BeamCalc::Instance()->getBandID(elfreq, telescopeName_p, bandName_p);
278 : }
279 :
280 0 : return bandID;
281 : };
282 0 : int EVLAAperture::getVisParams(const VisBuffer2 &vb,
283 : const CoordinateSystem & /*im*/) {
284 0 : throw(AipsError("EVLAAperture::getVisParams() called"));
285 : Double Freq;
286 : cacheVBInfo(vb);
287 :
288 : Freq = vb.getFrequency(0, 0);
289 : Double Lambda = C::c / Freq;
290 : HPBW = Lambda / (Diameter_p * sqrt(log(2.0)));
291 : sigma = 1.0 / (HPBW * HPBW);
292 :
293 : return getBandID(Freq);
294 : // Int bandID=0;
295 : // if (!isNoOp())
296 : // bandID = BeamCalc::Instance()->getBandID(Freq,telescopeName_p);
297 :
298 : // return bandID;
299 : }
300 :
301 0 : Int EVLAAperture::makePBPolnCoords(const VisBuffer2 &vb, const Int &convSize,
302 : const Int &convSampling,
303 : const CoordinateSystem &skyCoord,
304 : const Int &skyNx, const Int & /*skyNy*/,
305 : CoordinateSystem &feedCoord)
306 : // Vector<Int>& cfStokes)
307 : {
308 0 : feedCoord = skyCoord;
309 : //
310 : // Make a two dimensional image to calculate auto-correlation of
311 : // the ideal illumination pattern. We want this on a fine grid in
312 : // the UV plane
313 : //
314 0 : Int directionIndex = skyCoord.findCoordinate(Coordinate::DIRECTION);
315 : // cout<<"Direction index is "<< directionIndex<<endl;
316 0 : AlwaysAssert(directionIndex >= 0, AipsError);
317 0 : DirectionCoordinate dc = skyCoord.directionCoordinate(directionIndex);
318 0 : Vector<Double> sampling;
319 0 : sampling = dc.increment();
320 : // cout<<"Image sampling set to : "<<sampling<<endl;
321 0 : sampling *= Double(convSampling);
322 0 : sampling *= Double(skyNx) / Double(convSize);
323 0 : dc.setIncrement(sampling);
324 : // cout<<"Resized sampling is : "<<sampling<<endl;
325 :
326 0 : Vector<Double> unitVec(2);
327 0 : unitVec = convSize / 2;
328 0 : dc.setReferencePixel(unitVec);
329 :
330 : // Set the reference value to that of the image
331 0 : feedCoord.replaceCoordinate(dc, directionIndex);
332 :
333 : //
334 : // Make an image with circular polarization axis.
335 : //
336 0 : Int NPol = 0, M, N = 0;
337 0 : M = polMap_p_base.nelements();
338 0 : for (Int i = 0; i < M; i++)
339 0 : if (polMap_p_base(i) > -1)
340 0 : NPol++;
341 0 : Vector<Int> poln(NPol);
342 :
343 : Int index;
344 0 : Vector<Int> inStokes;
345 0 : index = feedCoord.findCoordinate(Coordinate::STOKES);
346 0 : inStokes = feedCoord.stokesCoordinate(index).stokes();
347 0 : N = 0;
348 : try {
349 : // cerr << "### " << polMap_p_base << " " << vb.corrType() << endl;
350 0 : for (Int i = 0; i < M; i++)
351 0 : if (polMap_p_base(i) > -1) {
352 0 : poln(N) = vb.correlationTypes()(i);
353 0 : N++;
354 : }
355 0 : StokesCoordinate polnCoord(poln);
356 0 : Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
357 0 : feedCoord.replaceCoordinate(polnCoord, StokesIndex);
358 : // cfStokes = poln;
359 0 : } catch (AipsError &x) {
360 0 : throw(SynthesisFTMachineError(
361 : "Likely cause: Discrepancy between the poln. "
362 0 : "axis of the data and the image specifications."));
363 0 : }
364 :
365 0 : return NPol;
366 0 : }
367 :
368 0 : Bool EVLAAperture::findSupport(Array<Complex> &func, Float &threshold,
369 : Int &origin, Int &R) {
370 : Double NSteps;
371 0 : Int PixInc = 1;
372 0 : Vector<Complex> vals;
373 0 : IPosition ndx(4, origin, 0, 0, 0);
374 0 : Bool found = false;
375 0 : IPosition cfShape = func.shape();
376 0 : Int convSize = cfShape(0);
377 0 : for (R = convSize / 4; R > 1; R--) {
378 0 : NSteps = 90 * R / PixInc; // Check every PixInc pixel along a
379 : // circle of radious R
380 0 : vals.resize((Int)(NSteps + 0.5));
381 0 : vals = 0;
382 0 : for (Int th = 0; th < NSteps; th++) {
383 0 : ndx(0) = (int)(origin + R * sin(2.0 * M_PI * th * PixInc / R));
384 0 : ndx(1) = (int)(origin + R * cos(2.0 * M_PI * th * PixInc / R));
385 :
386 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
387 0 : vals(th) = func(ndx);
388 : }
389 0 : if (max(abs(vals)) > threshold) {
390 0 : found = true;
391 0 : break;
392 : }
393 : }
394 0 : return found;
395 0 : }
396 :
397 0 : void EVLAAperture::makeFullJones(ImageInterface<Complex> &pbImage,
398 : const VisBuffer2 &vb, Bool doSquint,
399 : Int &bandID, Double freqVal) {
400 :
401 0 : if (!isNoOp()) {
402 0 : VLACalcIlluminationConvFunc vlaPB;
403 0 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
404 0 : vlaPB.setMaximumCacheSize(cachesize);
405 0 : bandID = getBandID(freqVal);
406 0 : vlaPB.makeFullJones(pbImage, vb, doSquint, bandID, freqVal);
407 0 : }
408 0 : }
409 :
410 0 : void EVLAAperture::applyAvgSkyJones(ImageInterface<Complex> &outImage) {
411 0 : TempImage<Complex> temp(outImage.shape(), outImage.coordinates());
412 0 : temp.setMiscInfo(outImage.miscInfo());
413 0 : temp.set(1.0);
414 0 : applyDiagSkyJones(temp, 0.0);
415 : // Taking the RL beam
416 0 : IPosition blc(4, 0, 0, 1, 0);
417 0 : IPosition trc=outImage.shape()-1;
418 0 : trc(2)=1;
419 0 : Slicer sl(blc, trc, Slicer::endIsLast);
420 0 : SubImage<Complex> rl(temp, sl, false);
421 0 : for (Int k = 0; k <4; ++k ) {
422 0 : trc[2] = k;
423 0 : blc[2] = k;
424 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
425 0 : SubImage<Complex> outsub(outImage, sl, true);
426 0 : outsub.copyData(LatticeExpr<Complex>(real(rl)));
427 0 : }
428 :
429 :
430 0 : }
431 :
432 0 : void EVLAAperture::applyDiagSkyJones(ImageInterface<Complex> &outImage,
433 : const Double pa) {
434 :
435 0 : ApertureCalcParams ap;
436 0 : Long memtot=HostInfo::memoryFree();
437 0 : Double memtobeused= Double(memtot)*1024.0;
438 0 : String bandname="";
439 : //cerr << "diagMisc " << outImage.miscInfo() << endl;
440 0 : if(outImage.miscInfo().isDefined("bandname"))
441 0 : outImage.miscInfo().get("bandname", bandname);
442 : //cerr << "diagSky BANDNAME " << bandname << endl;
443 :
444 0 : ap.apertureptr = make_shared< TempImage<Complex> >(outImage.shape(), outImage.coordinates(), memtobeused/10.0);
445 0 : ap.aperture = ap.apertureptr.get();
446 0 : CoordinateSystem csys = outImage.coordinates();
447 :
448 0 : Int index = csys.findCoordinate(Coordinate::SPECTRAL);
449 0 : Double freqVal = 0.0;
450 0 : csys.spectralCoordinate(index).toWorld(freqVal, 0.0);
451 : //cerr << "####FREQVAL " << freqVal << endl;
452 0 : Int bandID = getBandID(freqVal, bandname);
453 : //cerr << "bandid " << bandID << endl;
454 : CoordinateSystem uvCoords = refim::VLACalcIlluminationConvFunc::makeUVCoords(
455 0 : csys, outImage.shape(), freqVal);
456 0 : index = uvCoords.findCoordinate(Coordinate::LINEAR);
457 0 : if (index < 0)
458 0 : throw(AipsError("Could not get uv-coordinates in applyDiagSkyJones"));
459 0 : LinearCoordinate lc = uvCoords.linearCoordinate(index);
460 0 : Vector<Double> uvIncr = lc.increment();
461 0 : setApertureParams(ap, freqVal, pa, bandID, outImage.shape(), uvIncr);
462 : //ap.aperture->setCoordinateInfo(uvCoords);
463 : // Can replace this call with PBMathInterface::applyVP for other type of known beams
464 0 : BeamCalc::Instance()->calculateAperture(&ap);
465 : //cerr << "MAx MIN " << max(ap.aperture->get()) << " " << min(ap.aperture->get()) << endl;
466 : // need to ft to voltagepattern
467 0 : FFT2D ft;
468 0 : ft.c2cFFT(*(ap.apertureptr));
469 : // diagonal Mueller beam terms
470 : // First lets normalize voltage beam
471 : // In VLACalcIlluminationConvFunc::skyMuller function there is a double conj
472 : // first when normalizing tmp ...then when calculation M0 etc
473 : // So here we are just doing the conj once at multiplication
474 : // Copying the normalization from skyMuller ...hey it was tested in python must be right
475 0 : IPosition mid = outImage.shape()-1;
476 0 : mid[0] = outImage.shape()[0]/2;
477 0 : mid[1] = outImage.shape()[1]/2;
478 0 : Float normalizesq = 0.0;
479 0 : for (Int k = 0; k < 4; ++k ) {
480 0 : mid[2] = k;
481 0 : Complex valmid = ap.aperture->getAt(mid);
482 : //cerr << " mid " << mid << " val " << valmid << endl;
483 0 : normalizesq += abs(valmid*valmid)/2.0;
484 : }
485 0 : if (normalizesq == 0.0) {
486 0 : throw(AipsError("Voltage patterns are 0 at the center"));
487 : }
488 0 : (ap.aperture)-> copyData(LatticeExpr<Complex>((*(ap.aperture))/sqrt(normalizesq)));
489 : // now lets calculate the diag Mueller terms remember we are conj the multipliers to match the
490 : // double conj that VLACalcIlluminationConvFunc does !
491 : // For heterogenous arrays that is where we'll multiply the voltage pattern of ant1 with that
492 : // of ant2
493 :
494 : {
495 : // We'll overwrite the 2 middle planes with Jp*conj(Jq) and Jq*conj(Jp) (see skyMuller)
496 : // because we did not do a conj while normalizing it will be conj(Jp)*Jq and conj(Jq)*Jp
497 : // respectively
498 0 : IPosition blc(4, 0, 0, 0, 0);
499 0 : IPosition trc=outImage.shape()-1;
500 0 : trc(2)=0;
501 0 : Slicer sl(blc, trc, Slicer::endIsLast);
502 0 : SubImage<Complex> jp(*(ap.aperture), sl, true);
503 0 : blc(2) = 3;
504 0 : trc(2) = 3;
505 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
506 0 : SubImage<Complex> jq(*(ap.aperture), sl, true);
507 0 : blc(2) = 1;
508 0 : trc(2) = 1;
509 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
510 0 : SubImage<Complex> jpq(*(ap.aperture), sl, true);
511 0 : jpq.copyData(LatticeExpr<Complex>(jq*conj(jp)));
512 0 : blc(2) = 2;
513 0 : trc(2) = 2;
514 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
515 0 : SubImage<Complex> jqp(*(ap.aperture), sl, true);
516 0 : jqp.copyData(LatticeExpr<Complex>(conj(jq)*jp));
517 : // Now the first and last plane
518 0 : jp.copyData(LatticeExpr<Complex>(conj(jp)*jp));
519 0 : jq.copyData(LatticeExpr<Complex>(conj(jq)*jq));
520 0 : }
521 : //cerr << "SHAPES " << ap.aperture->shape() << " " << outImage.shape() << endl;
522 0 : LatticeExpr<Complex> le((*ap.aperture) * outImage);
523 0 : outImage.copyData(le);
524 0 : }
525 :
526 0 : void EVLAAperture::applySky(ImageInterface<Complex> &outImages,
527 : const VisBuffer2 &vb, const Bool doSquint,
528 : const Int &cfKey, const Int &muellerTerm,
529 : const Double freqVal) {
530 0 : Double freq = freqVal;
531 0 : if (freq < 0.0)
532 0 : freq = vb.getFrequency(0, 0);
533 : // (void)cfKey;
534 : // if (!isNoOp())
535 : // {
536 : // VLACalcIlluminationConvFunc vlaPB;
537 : // Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
538 : // vlaPB.setMaximumCacheSize(cachesize);
539 : // bandID = getBandID(freqVal,telescopeName_p);
540 : // // cout<<"EVLAAperture : muellerTerm"<<muellerTerm <<"\n";
541 : // //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
542 : // Double pa=getPA(vb);
543 : // vlaPB.applyPB(outImages, pa, doSquint,bandID,muellerTerm,freqVal);
544 :
545 : // }
546 0 : Double pa = getPA(vb);
547 : //cerr << "PA=" << pa << endl;
548 0 : applySky(outImages, pa, doSquint, cfKey, muellerTerm, freq);
549 0 : }
550 : //
551 : // This one without VB. Should become the default (and the only one!).
552 : //
553 0 : void EVLAAperture::applySky(ImageInterface<Complex> &outImages,
554 : const Double &pa, const Bool doSquint,
555 : const Int &cfKey, const Int &muellerTerm,
556 : const Double freqVal) {
557 : (void)cfKey;
558 0 : if (!isNoOp()) {
559 0 : VLACalcIlluminationConvFunc vlaPB;
560 0 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
561 0 : vlaPB.setMaximumCacheSize(cachesize);
562 : Int bandID;
563 0 : bandID = getBandID(freqVal);
564 : // vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
565 0 : Double pa_l = pa; // Due to goofup in making sure complier type checking
566 : // does not come in the way!
567 : //cerr << "PA_L=" << pa_l << " bandID= " << bandID << endl;
568 0 : vlaPB.applyPB(outImages, pa_l, doSquint, bandID, muellerTerm, freqVal);
569 0 : }
570 0 : }
571 :
572 0 : void EVLAAperture::applySky(ImageInterface<Float> &outImages,
573 : const VisBuffer2 &vb, const Bool doSquint,
574 : const Int &cfKey, const Int &muellerTerm,
575 : const Double freqVal) {
576 : (void)cfKey;
577 : (void)muellerTerm;
578 0 : Double freq = freqVal;
579 0 : if (freq < 0.0)
580 0 : freq = vb.getFrequency(0, 0);
581 0 : if (!isNoOp()) {
582 0 : VLACalcIlluminationConvFunc vlaPB;
583 0 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
584 0 : vlaPB.setMaximumCacheSize(cachesize);
585 : Int bandID;
586 0 : bandID = getBandID(freq);
587 :
588 0 : Double pa = getPA(vb);
589 0 : vlaPB.applyPB(outImages, pa, bandID, doSquint, freq);
590 0 : }
591 0 : }
592 :
593 : }; // namespace casa
|
