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 134 : EVLAAperture::EVLAAperture() : AzElAperture(), polMap_p(), feedStokes_p() {
55 134 : telescopeName_p = "EVLA";
56 134 : Diameter_p = 25.0;
57 134 : }
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 55 : void EVLAAperture::cacheVBInfo(const String &telescopeName,
155 : const Float &diameter) {
156 55 : telescopeName_p = telescopeName;
157 55 : Diameter_p = diameter;
158 55 : }
159 :
160 13 : void EVLAAperture::cacheVBInfo(const VisBuffer2 &vb) {
161 : const Vector<String> telescopeNames =
162 13 : vb.subtableColumns().observation().telescopeName().getColumn();
163 26 : for (uInt nt = 0; nt < telescopeNames.nelements(); nt++) {
164 13 : 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 13 : 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 13 : telescopeName_p = telescopeNames[0];
179 :
180 : // MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
181 : // Freq = vb.msColumns().spectralWindow().refFrequency()(0);
182 13 : Diameter_p = 0;
183 13 : Nant_p = vb.subtableColumns().antenna().nrow();
184 13 : for (Int i = 0; i < Nant_p; i++)
185 13 : if (!vb.subtableColumns().antenna().flagRow()(i)) {
186 13 : Diameter_p = vb.subtableColumns().antenna().dishDiameter().getColumn()(i);
187 13 : break;
188 : }
189 13 : 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 13 : cacheVBInfo(telescopeNames[0], Diameter_p);
196 13 : }
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 : bandID = BeamCalc::Instance()->getBandID(freq, telescopeName_p, tokens(0));
206 0 : }
207 :
208 0 : return bandID;
209 : };
210 148 : String EVLAAperture::getVLABandName(const Double& freq, const String& telescopeName) {
211 :
212 148 : String bandName = "EVLA_L";
213 148 : if (telescopeName == "VLA") {
214 : // if ((freq >= 1.34E9) && (freq <= 1.73E9))
215 0 : if((freq >= 9E8) && (freq <= 1.73E9))
216 0 : bandName = "VLA_L";
217 0 : else if ((freq >= 4.5E9) && (freq <= 5.0E9))
218 0 : bandName = "VLA_C";
219 0 : else if ((freq >= 8.0E9) && (freq <= 8.8E9))
220 0 : bandName = "VLA_X";
221 0 : else if ((freq >= 14.4E9) && (freq <= 15.4E9))
222 0 : bandName = "VLA_U";
223 0 : else if ((freq >= 22.0E9) && (freq <= 24.0E9))
224 0 : bandName = "VLA_K";
225 0 : else if ((freq >= 40.0E9) && (freq <= 50.0E9))
226 0 : bandName = "VLA_Q";
227 0 : else if ((freq >= 30E6) && (freq <= 100E6))
228 0 : bandName = "VLA_4";
229 : else
230 0 : throw(
231 0 : AipsError("Don't know VLA band for frequency=" + String::toString(freq)));
232 148 : } else if (telescopeName == "EVLA") {
233 148 : if (freq > 9e8 && freq <= 2.0e9)
234 92 : bandName = "EVLA_L";
235 56 : else if (freq > 2.0e9 && freq <= 4.0e9)
236 56 : bandName = "EVLA_S";
237 0 : else if (freq > 4.0e9 && freq <= 8.0e9)
238 0 : bandName = "EVLA_C";
239 0 : else if (freq > 8.e9 && freq <= 12.0e9)
240 0 : bandName = "EVLA_X";
241 0 : else if (freq > 12.0e9 && freq <= 18.0e9)
242 0 : bandName = "EVLA_U";
243 0 : else if (freq > 18.0e9 && freq <= 26.0e9)
244 0 : bandName = "EVLA_K";
245 0 : else if (freq > 26.e9 && freq <= 40.0e9)
246 0 : bandName = "EVLA_A";
247 0 : else if (freq > 40.0e9 && freq <= 50.0e9)
248 0 : bandName = "EVLA_Q";
249 : else
250 0 : throw(
251 0 : AipsError("Don't know EVLA band for frequency=" + String::toString(freq)));
252 : } else {
253 0 : throw(AipsError("Don't know telescope " + telescopeName));
254 : }
255 148 : return bandName;
256 0 : }
257 148 : Int EVLAAperture::getBandID(const Double &freq, const String& bandName) {
258 148 : Int bandID = 0;
259 148 : if(bandName=="")
260 148 : bandName_p = getVLABandName(freq, telescopeName_p);
261 : else
262 0 : bandName_p=bandName;
263 :
264 148 : if (!isNoOp()) {
265 : // First #-separated token in bandName_p is the name of the band used
266 : // Vector<String> tokens = SynthesisUtils::parseBandName(bandName_p);
267 : // cerr << "TOKENS " << tokens << endl;
268 148 : bandID = BeamCalc::Instance()->getBandID(freq, telescopeName_p, bandName_p);
269 : }
270 :
271 148 : return bandID;
272 : };
273 0 : int EVLAAperture::getVisParams(const VisBuffer2 &vb,
274 : const CoordinateSystem & /*im*/) {
275 0 : throw(AipsError("EVLAAperture::getVisParams() called"));
276 : Double Freq;
277 : cacheVBInfo(vb);
278 :
279 : Freq = vb.getFrequency(0, 0);
280 : Double Lambda = C::c / Freq;
281 : HPBW = Lambda / (Diameter_p * sqrt(log(2.0)));
282 : sigma = 1.0 / (HPBW * HPBW);
283 :
284 : return getBandID(Freq);
285 : // Int bandID=0;
286 : // if (!isNoOp())
287 : // bandID = BeamCalc::Instance()->getBandID(Freq,telescopeName_p);
288 :
289 : // return bandID;
290 : }
291 :
292 13 : Int EVLAAperture::makePBPolnCoords(const VisBuffer2 &vb, const Int &convSize,
293 : const Int &convSampling,
294 : const CoordinateSystem &skyCoord,
295 : const Int &skyNx, const Int & /*skyNy*/,
296 : CoordinateSystem &feedCoord)
297 : // Vector<Int>& cfStokes)
298 : {
299 13 : feedCoord = skyCoord;
300 : //
301 : // Make a two dimensional image to calculate auto-correlation of
302 : // the ideal illumination pattern. We want this on a fine grid in
303 : // the UV plane
304 : //
305 13 : Int directionIndex = skyCoord.findCoordinate(Coordinate::DIRECTION);
306 : // cout<<"Direction index is "<< directionIndex<<endl;
307 13 : AlwaysAssert(directionIndex >= 0, AipsError);
308 13 : DirectionCoordinate dc = skyCoord.directionCoordinate(directionIndex);
309 13 : Vector<Double> sampling;
310 13 : sampling = dc.increment();
311 : // cout<<"Image sampling set to : "<<sampling<<endl;
312 13 : sampling *= Double(convSampling);
313 13 : sampling *= Double(skyNx) / Double(convSize);
314 13 : dc.setIncrement(sampling);
315 : // cout<<"Resized sampling is : "<<sampling<<endl;
316 :
317 13 : Vector<Double> unitVec(2);
318 13 : unitVec = convSize / 2;
319 13 : dc.setReferencePixel(unitVec);
320 :
321 : // Set the reference value to that of the image
322 13 : feedCoord.replaceCoordinate(dc, directionIndex);
323 :
324 : //
325 : // Make an image with circular polarization axis.
326 : //
327 13 : Int NPol = 0, M, N = 0;
328 13 : M = polMap_p_base.nelements();
329 39 : for (Int i = 0; i < M; i++)
330 26 : if (polMap_p_base(i) > -1)
331 26 : NPol++;
332 13 : Vector<Int> poln(NPol);
333 :
334 : Int index;
335 13 : Vector<Int> inStokes;
336 13 : index = feedCoord.findCoordinate(Coordinate::STOKES);
337 13 : inStokes = feedCoord.stokesCoordinate(index).stokes();
338 13 : N = 0;
339 : try {
340 : // cerr << "### " << polMap_p_base << " " << vb.corrType() << endl;
341 39 : for (Int i = 0; i < M; i++)
342 26 : if (polMap_p_base(i) > -1) {
343 26 : poln(N) = vb.correlationTypes()(i);
344 26 : N++;
345 : }
346 13 : StokesCoordinate polnCoord(poln);
347 13 : Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
348 13 : feedCoord.replaceCoordinate(polnCoord, StokesIndex);
349 : // cfStokes = poln;
350 13 : } catch (AipsError &x) {
351 0 : throw(SynthesisFTMachineError(
352 : "Likely cause: Discrepancy between the poln. "
353 0 : "axis of the data and the image specifications."));
354 0 : }
355 :
356 13 : return NPol;
357 13 : }
358 :
359 0 : Bool EVLAAperture::findSupport(Array<Complex> &func, Float &threshold,
360 : Int &origin, Int &R) {
361 : Double NSteps;
362 0 : Int PixInc = 1;
363 0 : Vector<Complex> vals;
364 0 : IPosition ndx(4, origin, 0, 0, 0);
365 0 : Bool found = false;
366 0 : IPosition cfShape = func.shape();
367 0 : Int convSize = cfShape(0);
368 0 : for (R = convSize / 4; R > 1; R--) {
369 0 : NSteps = 90 * R / PixInc; // Check every PixInc pixel along a
370 : // circle of radious R
371 0 : vals.resize((Int)(NSteps + 0.5));
372 0 : vals = 0;
373 0 : for (Int th = 0; th < NSteps; th++) {
374 0 : ndx(0) = (int)(origin + R * sin(2.0 * M_PI * th * PixInc / R));
375 0 : ndx(1) = (int)(origin + R * cos(2.0 * M_PI * th * PixInc / R));
376 :
377 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
378 0 : vals(th) = func(ndx);
379 : }
380 0 : if (max(abs(vals)) > threshold) {
381 0 : found = true;
382 0 : break;
383 : }
384 : }
385 0 : return found;
386 0 : }
387 :
388 0 : void EVLAAperture::makeFullJones(ImageInterface<Complex> &pbImage,
389 : const VisBuffer2 &vb, Bool doSquint,
390 : Int &bandID, Double freqVal) {
391 :
392 0 : if (!isNoOp()) {
393 0 : VLACalcIlluminationConvFunc vlaPB;
394 0 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
395 0 : vlaPB.setMaximumCacheSize(cachesize);
396 0 : bandID = getBandID(freqVal);
397 0 : vlaPB.makeFullJones(pbImage, vb, doSquint, bandID, freqVal);
398 0 : }
399 0 : }
400 :
401 0 : void EVLAAperture::applyAvgSkyJones(ImageInterface<Complex> &outImage) {
402 0 : TempImage<Complex> temp(outImage.shape(), outImage.coordinates());
403 0 : temp.setMiscInfo(outImage.miscInfo());
404 0 : temp.set(1.0);
405 0 : applyDiagSkyJones(temp, 0.0);
406 : // Taking the RL beam
407 0 : IPosition blc(4, 0, 0, 1, 0);
408 0 : IPosition trc=outImage.shape()-1;
409 0 : trc(2)=1;
410 0 : Slicer sl(blc, trc, Slicer::endIsLast);
411 0 : SubImage<Complex> rl(temp, sl, false);
412 0 : for (Int k = 0; k <4; ++k ) {
413 0 : trc[2] = k;
414 0 : blc[2] = k;
415 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
416 0 : SubImage<Complex> outsub(outImage, sl, true);
417 0 : outsub.copyData(LatticeExpr<Complex>(real(rl)));
418 0 : }
419 :
420 :
421 0 : }
422 :
423 0 : void EVLAAperture::applyDiagSkyJones(ImageInterface<Complex> &outImage,
424 : const Double pa) {
425 :
426 0 : ApertureCalcParams ap;
427 0 : Long memtot=HostInfo::memoryFree();
428 0 : Double memtobeused= Double(memtot)*1024.0;
429 0 : String bandname="";
430 : //cerr << "diagMisc " << outImage.miscInfo() << endl;
431 0 : if(outImage.miscInfo().isDefined("bandname"))
432 0 : outImage.miscInfo().get("bandname", bandname);
433 : //cerr << "diagSky BANDNAME " << bandname << endl;
434 :
435 0 : ap.apertureptr = make_shared< TempImage<Complex> >(outImage.shape(), outImage.coordinates(), memtobeused/10.0);
436 0 : ap.aperture = ap.apertureptr.get();
437 0 : CoordinateSystem csys = outImage.coordinates();
438 :
439 0 : Int index = csys.findCoordinate(Coordinate::SPECTRAL);
440 0 : Double freqVal = 0.0;
441 0 : csys.spectralCoordinate(index).toWorld(freqVal, 0.0);
442 : //cerr << "####FREQVAL " << freqVal << endl;
443 0 : Int bandID = getBandID(freqVal, bandname);
444 : //cerr << "bandid " << bandID << endl;
445 : CoordinateSystem uvCoords = refim::VLACalcIlluminationConvFunc::makeUVCoords(
446 0 : csys, outImage.shape(), freqVal);
447 0 : index = uvCoords.findCoordinate(Coordinate::LINEAR);
448 0 : if (index < 0)
449 0 : throw(AipsError("Could not get uv-coordinates in applyDiagSkyJones"));
450 0 : LinearCoordinate lc = uvCoords.linearCoordinate(index);
451 0 : Vector<Double> uvIncr = lc.increment();
452 0 : setApertureParams(ap, freqVal, pa, bandID, outImage.shape(), uvIncr);
453 : //ap.aperture->setCoordinateInfo(uvCoords);
454 : // Can replace this call with PBMathInterface::applyVP for other type of known beams
455 0 : BeamCalc::Instance()->calculateAperture(&ap);
456 : //cerr << "MAx MIN " << max(ap.aperture->get()) << " " << min(ap.aperture->get()) << endl;
457 : // need to ft to voltagepattern
458 0 : FFT2D ft;
459 0 : ft.c2cFFT(*(ap.apertureptr));
460 : // diagonal Mueller beam terms
461 : // First lets normalize voltage beam
462 : // In VLACalcIlluminationConvFunc::skyMuller function there is a double conj
463 : // first when normalizing tmp ...then when calculation M0 etc
464 : // So here we are just doing the conj once at multiplication
465 : // Copying the normalization from skyMuller ...hey it was tested in python must be right
466 0 : IPosition mid = outImage.shape()-1;
467 0 : mid[0] = outImage.shape()[0]/2;
468 0 : mid[1] = outImage.shape()[1]/2;
469 0 : Float normalizesq = 0.0;
470 0 : for (Int k = 0; k < 4; ++k ) {
471 0 : mid[2] = k;
472 0 : Complex valmid = ap.aperture->getAt(mid);
473 : //cerr << " mid " << mid << " val " << valmid << endl;
474 0 : normalizesq += abs(valmid*valmid)/2.0;
475 : }
476 0 : if (normalizesq == 0.0) {
477 0 : throw(AipsError("Voltage patterns are 0 at the center"));
478 : }
479 0 : (ap.aperture)-> copyData(LatticeExpr<Complex>((*(ap.aperture))/sqrt(normalizesq)));
480 : // now lets calculate the diag Mueller terms remember we are conj the multipliers to match the
481 : // double conj that VLACalcIlluminationConvFunc does !
482 : // For heterogenous arrays that is where we'll multiply the voltage pattern of ant1 with that
483 : // of ant2
484 :
485 : {
486 : // We'll overwrite the 2 middle planes with Jp*conj(Jq) and Jq*conj(Jp) (see skyMuller)
487 : // because we did not do a conj while normalizing it will be conj(Jp)*Jq and conj(Jq)*Jp
488 : // respectively
489 0 : IPosition blc(4, 0, 0, 0, 0);
490 0 : IPosition trc=outImage.shape()-1;
491 0 : trc(2)=0;
492 0 : Slicer sl(blc, trc, Slicer::endIsLast);
493 0 : SubImage<Complex> jp(*(ap.aperture), sl, true);
494 0 : blc(2) = 3;
495 0 : trc(2) = 3;
496 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
497 0 : SubImage<Complex> jq(*(ap.aperture), sl, true);
498 0 : blc(2) = 1;
499 0 : trc(2) = 1;
500 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
501 0 : SubImage<Complex> jpq(*(ap.aperture), sl, true);
502 0 : jpq.copyData(LatticeExpr<Complex>(jq*conj(jp)));
503 0 : blc(2) = 2;
504 0 : trc(2) = 2;
505 0 : sl = Slicer(blc, trc, Slicer::endIsLast);
506 0 : SubImage<Complex> jqp(*(ap.aperture), sl, true);
507 0 : jqp.copyData(LatticeExpr<Complex>(conj(jq)*jp));
508 : // Now the first and last plane
509 0 : jp.copyData(LatticeExpr<Complex>(conj(jp)*jp));
510 0 : jq.copyData(LatticeExpr<Complex>(conj(jq)*jq));
511 0 : }
512 : //cerr << "SHAPES " << ap.aperture->shape() << " " << outImage.shape() << endl;
513 0 : LatticeExpr<Complex> le((*ap.aperture) * outImage);
514 0 : outImage.copyData(le);
515 0 : }
516 :
517 64 : void EVLAAperture::applySky(ImageInterface<Complex> &outImages,
518 : const VisBuffer2 &vb, const Bool doSquint,
519 : const Int &cfKey, const Int &muellerTerm,
520 : const Double freqVal) {
521 64 : Double freq = freqVal;
522 64 : if (freq < 0.0)
523 0 : freq = vb.getFrequency(0, 0);
524 : // (void)cfKey;
525 : // if (!isNoOp())
526 : // {
527 : // VLACalcIlluminationConvFunc vlaPB;
528 : // Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
529 : // vlaPB.setMaximumCacheSize(cachesize);
530 : // bandID = getBandID(freqVal,telescopeName_p);
531 : // // cout<<"EVLAAperture : muellerTerm"<<muellerTerm <<"\n";
532 : // //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
533 : // Double pa=getPA(vb);
534 : // vlaPB.applyPB(outImages, pa, doSquint,bandID,muellerTerm,freqVal);
535 :
536 : // }
537 64 : Double pa = getPA(vb);
538 : //cerr << "PA=" << pa << endl;
539 64 : applySky(outImages, pa, doSquint, cfKey, muellerTerm, freq);
540 64 : }
541 : //
542 : // This one without VB. Should become the default (and the only one!).
543 : //
544 148 : void EVLAAperture::applySky(ImageInterface<Complex> &outImages,
545 : const Double &pa, const Bool doSquint,
546 : const Int &cfKey, const Int &muellerTerm,
547 : const Double freqVal) {
548 : (void)cfKey;
549 148 : if (!isNoOp()) {
550 148 : VLACalcIlluminationConvFunc vlaPB;
551 148 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
552 148 : vlaPB.setMaximumCacheSize(cachesize);
553 : Int bandID;
554 148 : bandID = getBandID(freqVal);
555 : // vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
556 148 : Double pa_l = pa; // Due to goofup in making sure complier type checking
557 : // does not come in the way!
558 : //cerr << "PA_L=" << pa_l << " bandID= " << bandID << endl;
559 148 : vlaPB.applyPB(outImages, pa_l, doSquint, bandID, muellerTerm, freqVal);
560 148 : }
561 148 : }
562 :
563 0 : void EVLAAperture::applySky(ImageInterface<Float> &outImages,
564 : const VisBuffer2 &vb, const Bool doSquint,
565 : const Int &cfKey, const Int &muellerTerm,
566 : const Double freqVal) {
567 : (void)cfKey;
568 : (void)muellerTerm;
569 0 : Double freq = freqVal;
570 0 : if (freq < 0.0)
571 0 : freq = vb.getFrequency(0, 0);
572 0 : if (!isNoOp()) {
573 0 : VLACalcIlluminationConvFunc vlaPB;
574 0 : Long cachesize = (HostInfo::memoryTotal(true) / 8) * 1024;
575 0 : vlaPB.setMaximumCacheSize(cachesize);
576 : Int bandID;
577 0 : bandID = getBandID(freq);
578 :
579 0 : Double pa = getPA(vb);
580 0 : vlaPB.applyPB(outImages, pa, bandID, doSquint, freq);
581 0 : }
582 0 : }
583 :
584 : }; // namespace casa
|
