Line data Source code
1 : //# ALMAAperture.cc: Implementation of the ALMAAperture class
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //# Copyright (C) 2011 by ESO (in the framework of the ALMA collaboration)
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 <synthesis/TransformMachines/Utils.h>
30 : #include <synthesis/TransformMachines/ALMAAperture.h>
31 : #include <synthesis/TransformMachines/SynthesisError.h>
32 : #include <synthesis/TransformMachines/BeamCalc.h>
33 : #include <synthesis/TransformMachines/WTerm.h>
34 : #include <synthesis/TransformMachines/ALMACalcIlluminationConvFunc.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 :
39 : using namespace casacore;
40 : namespace casa{
41 :
42 : AntennaResponses* ALMAAperture::aR_p = 0;
43 : Bool ALMAAperture::orderMattersInCFKey = false;
44 :
45 0 : ALMAAperture::ALMAAperture():
46 : // ATerm(),
47 : AzElAperture(),
48 0 : polMap_p(),
49 0 : antTypeMap_p(),
50 0 : respImage_p()
51 : {
52 :
53 0 : haveCannedResponses_p = true;
54 :
55 0 : if(!aR_p){
56 0 : aR_p = new AntennaResponses();
57 : }
58 :
59 0 : if(!aR_p->isInit()){ // the shared antenna responses object was not yet initialised
60 :
61 0 : cout << "Initialising antenna responses ..." << endl;
62 :
63 0 : String antRespPath;
64 0 : if(!MeasTable::AntennaResponsesPath(antRespPath, "ALMA")) {
65 : // unknown observatory
66 0 : logIO() << LogOrigin("ALMAAperture", "ctor")
67 : << LogIO::WARN
68 : << "We don't have any precalculated antenna responses for ALMA."
69 : << endl << "Will try to generate them using ray tracing instead."
70 0 : << LogIO::POST;
71 0 : haveCannedResponses_p = false;
72 : }
73 0 : else if(!aR_p->init(antRespPath)){
74 : // init failed
75 0 : String mesg="Initialisation of antenna responses for observatory ALMA failed using path "+antRespPath;
76 0 : SynthesisError err(mesg);
77 0 : throw(err);
78 0 : }
79 0 : else if(MeasTable::AntennaResponsesPath(antRespPath, "ACA")) {
80 : // also have responses for ACA
81 0 : aR_p->append(antRespPath); // dont't throw if this fails because the ACA antennas
82 : // may already be in the ALMA table
83 : }
84 0 : }
85 :
86 0 : }
87 :
88 0 : ALMAAperture::~ALMAAperture(){
89 0 : for(uInt i=0;i<respImage_p.nelements();i++){
90 0 : if(respImage_p(i)){
91 0 : delete respImage_p(i);
92 : }
93 : }
94 0 : };
95 :
96 :
97 0 : ALMAAperture& ALMAAperture::operator=(const ALMAAperture& other)
98 : {
99 0 : if(this!=&other)
100 : {
101 : // aR_p does not need to be copied because it is static
102 0 : haveCannedResponses_p = other.haveCannedResponses_p;
103 0 : polMap_p.assign(other.polMap_p);
104 0 : antTypeMap_p.assign(other.antTypeMap_p);
105 : }
106 0 : return *this;
107 : }
108 0 : void ALMAAperture::cacheVBInfo(const String& telescopeName,
109 : const Float& diameter)
110 : {
111 0 : telescopeName_p=telescopeName;
112 0 : Diameter_p=diameter;
113 0 : }
114 0 : void ALMAAperture::cacheVBInfo(const VisBuffer& vb)
115 : {
116 0 : Vector<String> telescopeNames=vb.msColumns().observation().telescopeName().getColumn();
117 0 : for(uInt nt=0;nt<telescopeNames.nelements();nt++){
118 0 : if ((telescopeNames(nt) != "ALMA") &&
119 0 : (telescopeNames(nt) != "ACA") &&
120 0 : (telescopeNames(nt) != "OSF")){
121 0 : String mesg="Can handle only ALMA, ACA, and OSF antennas.\n";
122 0 : mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
123 0 : SynthesisError err(mesg);
124 0 : throw(err);
125 0 : }
126 : }
127 0 : telescopeName_p=telescopeNames[0];
128 0 : }
129 :
130 0 : Int ALMAAperture::getBandID(const Double& freq, const String& /*telescopeName*/)
131 : {
132 0 : Int bandID = -1;
133 0 : if(haveCannedResponses_p){
134 0 : String bandName;
135 0 : if(aR_p->getBandName(bandName, "ALMA", freq)){
136 0 : bandID = atoi(bandName.substr(bandName.find("and")+3).c_str()); // band names start with "band"
137 : }
138 : else{
139 0 : logIO() << LogOrigin("ALMAAperture", "getVisParams")
140 : << LogIO::WARN
141 : << "We don't have predefined antenna responses for ALMA at "
142 : << freq << " Hz. Will try to use raytracing instead."
143 0 : << LogIO::POST;
144 : }
145 0 : }
146 0 : return bandID;
147 : }
148 :
149 0 : Int ALMAAperture::getVisParams(const VisBuffer& vb,
150 : const CoordinateSystem& /*skyCoord*/)
151 : {
152 0 : cacheVBInfo(vb);
153 :
154 0 : MVFrequency FreqQ(vb.msColumns().spectralWindow().refFrequencyQuant()(0));
155 0 : Double Freq = FreqQ.getValue();
156 :
157 0 : return getBandID(Freq, telescopeName_p);
158 :
159 : // Int bandID = -1;
160 : // if(haveCannedResponses_p){
161 : // String bandName;
162 : // if(aR_p->getBandName(bandName, "ALMA", FreqQ)){
163 : // bandID = atoi(bandName.substr(bandName.find("and")+3).c_str()); // band names start with "band"
164 : // }
165 : // else{
166 : // logIO() << LogOrigin("ALMAAperture", "getVisParams")
167 : // << LogIO::WARN
168 : // << "We don't have predefined antenna responses for ALMA at "
169 : // << Freq << " Hz. Will try to use raytracing instead."
170 : // << LogIO::POST;
171 : // }
172 : // }
173 : // return bandID;
174 0 : }
175 :
176 0 : void ALMAAperture::makeFullJones(ImageInterface<Complex>& /*pbImage*/,
177 : const VisBuffer& /*vb*/,
178 : Bool /*doSquint*/, Int& /*bandID*/, Double /*freqVal*/)
179 : {
180 0 : throw(AipsError("ALMAAperture::makeJones() not yet implemented"));
181 : }
182 :
183 0 : void ALMAAperture::applySky(ImageInterface<Complex>& outImage, // the image (cube if there is more than one pol)
184 : // upon return: FTed convolution pair of AIFs given by the cfKey
185 : // rotated by the PA and regridded to the image
186 : const VisBuffer& vb, // for the parallactic angle
187 : const Bool doSquint,
188 : const Int& cfKey,
189 : const Bool raytrace)
190 : {
191 :
192 0 : if(getVisParams(vb)==-1 || raytrace){ // need to use ray tracing
193 0 : ALMACalcIlluminationConvFunc almaPB;
194 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
195 0 : almaPB.setMaximumCacheSize(cachesize);
196 0 : almaPB.applyPB(outImage, vb, doSquint, cfKey);
197 0 : }
198 : else{ // use canned antenna responses
199 :
200 : // (issue warning if image too coarse => will be in method ATerm::OK(vb, PAtolerance, timetolerance)
201 :
202 : // identify the right response image based on antenna, freq band, polarizations
203 0 : Vector<ALMAAntennaType> aTypes = antennaTypesFromCFKey(cfKey);
204 0 : uInt nAntTypes = aTypes.nelements();
205 0 : Int spwId = vb.spectralWindow();
206 :
207 0 : Vector<String> respImageName(nAntTypes);
208 0 : Vector<uInt> respImageChannel(nAntTypes);
209 0 : Vector<MFrequency> respImageNomFreq(nAntTypes);
210 0 : Vector<AntennaResponses::FuncTypes> respImageFType(nAntTypes);
211 0 : Vector<MVAngle> respImageRotOffset(nAntTypes);
212 0 : Vector<Vector<Array<Complex> > > respByPol(nAntTypes);
213 :
214 0 : MFrequency refFreq = vb.msColumns().spectralWindow().refFrequencyMeas()(spwId);
215 0 : MEpoch obsTime = vb.msColumns().observation().timeRangeMeas()(vb.observationId()(0))(IPosition(1,0)); // start time
216 0 : MDirection obsDir = vb.msColumns().field().phaseDirMeas(vb.fieldId());
217 :
218 0 : for(uInt i=0; i<nAntTypes; i++){
219 0 : cout << "Looking for ant type \"" << antTypeStrFromType(aTypes(i)) << "\", freq "
220 0 : << refFreq.get(Unit("Hz")).getValue() << " Hz" << endl;
221 0 : if(!aR_p->getImageName(respImageName(i),
222 0 : respImageChannel(i),
223 : respImageNomFreq(i),
224 : respImageFType(i),
225 : respImageRotOffset(i),
226 : "ALMA",
227 : obsTime,
228 : refFreq,
229 0 : AntennaResponses::EFP,
230 0 : antTypeStrFromType(aTypes(i)),
231 : obsDir,
232 : "", // receiver type
233 0 : 0 // beam number
234 : )){ // no matching response found
235 0 : ostringstream oss;
236 0 : oss << "No matching antenna response found for frequency "
237 0 : << refFreq.get(Unit("Hz")).getValue() << " Hz, and antenna type "
238 0 : << antTypeStrFromType(aTypes(i));
239 0 : throw(SynthesisError(oss.str()));
240 0 : }
241 : }
242 :
243 : // load all necessary images (max. two)
244 0 : for(uInt i=0; i<nAntTypes; i++){
245 0 : if(respImageFType(i)!=AntennaResponses::EFP){
246 0 : throw(SynthesisError(String("Can only process antenna responses of type EFP.")));
247 : }
248 : }
249 0 : respImage_p.resize(nAntTypes);
250 0 : for(uInt i=0; i<nAntTypes; i++){
251 0 : cout << "Loading " << respImageName(i) << endl;
252 :
253 : try{
254 0 : respImage_p(i) = new PagedImage<Complex>(respImageName(i));
255 : }
256 0 : catch(std::exception x){
257 0 : ostringstream oss;
258 0 : oss << "Error reading antenna response image from path \""
259 0 : << respImageName(i) << "\": " << x.what();
260 0 : respImage_p.resize(i,true);
261 0 : throw(SynthesisError(oss.str()));
262 0 : }
263 : }
264 0 : cout << "Loaded " << nAntTypes << " images." << endl;
265 :
266 :
267 :
268 : // check if there are spectral and stokes axes in the output image
269 0 : Int pIndex = outImage.coordinates().findCoordinate(Coordinate::STOKES);
270 0 : if(pIndex==-1){
271 0 : ostringstream oss;
272 0 : oss << "Error: input image does not contain the necessary polarisation axis.";
273 0 : throw(SynthesisError(oss.str()));
274 0 : }
275 0 : Int pSIndex = outImage.coordinates().findCoordinate(Coordinate::SPECTRAL);
276 :
277 : // identify polarisations in the response image and fill into respByPol(i)
278 0 : Int rIndex = respImage_p(0)->coordinates().findCoordinate(Coordinate::STOKES);
279 0 : Int rAxis = respImage_p(0)->coordinates().pixelAxes(rIndex)(0);
280 0 : Vector<Int> rStokes = respImage_p(0)->coordinates().stokesCoordinate(rIndex).stokes();
281 0 : IPosition rShape = respImage_p(0)->shape();
282 0 : if( (nAntTypes>1) && rShape != respImage_p(1)->shape()){
283 0 : ostringstream oss;
284 0 : oss << "Error: response images for different antenna types (but otherwise identical parameters)"
285 0 : << endl << "need to have the same shape:"
286 0 : << "Resp. image 1 shape " << rShape << endl
287 0 : << "Resp. image 2 shape " << respImage_p(1)->shape() << endl;
288 0 : throw(SynthesisError(oss.str()));
289 0 : }
290 :
291 0 : const uInt rNDim = rShape.size();
292 :
293 0 : IPosition rSkyShape(rNDim,1);
294 0 : rSkyShape(0) = rShape(0);
295 0 : rSkyShape(1) = rShape(1);
296 0 : uInt nRPol = rStokes.nelements();
297 0 : if(!(nRPol==4 &&
298 0 : rStokes(0)==Stokes::XX && rStokes(1)==Stokes::XY && rStokes(2)==Stokes::YX && rStokes(3)==Stokes::YY)){
299 0 : ostringstream oss;
300 0 : oss << "Error: Antenna response image from path \""
301 0 : << respImageName(0)
302 : << "\" does not contain the necessary polarisation products or they are in the wrong order.\n"
303 0 : << "Order should be XX, XY, YX, YY.";
304 0 : throw(SynthesisError(oss.str()));
305 0 : }
306 0 : for(uInt i=0;i<nAntTypes; i++){
307 0 : respByPol(i).resize(nRPol);
308 : }
309 0 : for(uInt i=0; i<nAntTypes; i++){
310 0 : for(uInt iPol=0; iPol<nRPol; iPol++){
311 0 : IPosition start(rNDim,0);
312 0 : start(rAxis) = iPol;
313 0 : respImage_p(i)->getSlice(respByPol(i)(iPol), start, rSkyShape);
314 0 : }
315 : }
316 :
317 : // identify polarizations in the input image, put them into polToDoIndex()
318 0 : Vector<Int> inStokes;
319 0 : inStokes = outImage.coordinates().stokesCoordinate(pIndex).stokes();
320 0 : uInt nPol = inStokes.nelements();
321 0 : Vector<uInt> polToDoIndex(nPol);
322 0 : for(uInt i=0; i<nPol; i++){
323 0 : uInt ival=-1;
324 0 : switch(inStokes(i)){
325 0 : case Stokes::XX:
326 0 : ival = 0;
327 0 : break;
328 0 : case Stokes::XY:
329 0 : ival = 1;
330 0 : break;
331 0 : case Stokes::YX:
332 0 : ival = 2;
333 0 : break;
334 0 : case Stokes::YY:
335 0 : ival = 3;
336 0 : break;
337 0 : default:
338 0 : ostringstream oss;
339 0 : oss << "Error processing input image: polarization not valid for ALMA: "
340 0 : << Stokes::name(Stokes::type(inStokes(i)));
341 0 : throw(SynthesisError(oss.str()));
342 : break;
343 : }
344 0 : polToDoIndex(i) = ival;
345 : }
346 :
347 : // Calculate the primary beam for the given baseline type for each polarization
348 0 : LogIO os(LogOrigin("ALMAAperture", "applySky", WHERE));
349 :
350 0 : CoordinateSystem dCoord = respImage_p(0)->coordinates(); // assume both response images have same coordsys
351 :
352 0 : CoordinateSystem dCoordFinal(dCoord);
353 0 : uInt rNDimFinal = rNDim;
354 0 : IPosition dShapeFinal = respImage_p(0)->shape();
355 0 : dShapeFinal(rAxis) = nPol; // set the number of stokes pixels to that of the output image
356 :
357 : // check if we need to add a degenerate spectral axis to the response image
358 0 : Int rSIndex = dCoord.findCoordinate(Coordinate::SPECTRAL);
359 0 : if(rSIndex==-1 && pSIndex!=-1){// no spectral coordinate in resp. image but input has one
360 0 : SpectralCoordinate sC;
361 0 : dCoordFinal.addCoordinate(sC);
362 0 : rSIndex = dCoordFinal.findCoordinate(Coordinate::SPECTRAL);
363 0 : dShapeFinal.resize(4,true);
364 0 : dShapeFinal(3)=1;
365 0 : rNDimFinal +=1;
366 0 : }
367 0 : else if(rSIndex!=-1 && pSIndex==-1){ // no spectral coordinate in input image but response has one
368 0 : dCoordFinal = CoordinateSystem();
369 0 : dCoordFinal.addCoordinate(dCoord.directionCoordinate(dCoord.findCoordinate(Coordinate::DIRECTION)));
370 0 : dCoordFinal.addCoordinate(dCoord.stokesCoordinate(dCoord.findCoordinate(Coordinate::STOKES)));
371 0 : dShapeFinal.resize(3,true);
372 0 : rNDimFinal -=1;
373 0 : rSIndex = -1;
374 : }
375 :
376 0 : TempImage<Complex> nearFinal(dShapeFinal, dCoordFinal);
377 :
378 0 : for(uInt iPol=0; iPol<nPol; iPol++){
379 :
380 0 : Array<Complex> pB( respByPol(0)(polToDoIndex(iPol)).shape() );
381 0 : Array<Complex> fact1( pB.shape() );
382 0 : Array<Complex> fact2( pB.shape() );
383 :
384 : // rotate the two factor arrays into the right PA
385 0 : Double dAngleRad = getPA(vb);
386 0 : cout << "PA = " << dAngleRad << " rad" << endl;
387 :
388 0 : Int fact1Index=-1, fact2Index=-1;
389 : Double pA1, pA2;
390 :
391 : // apply the rotation offset from the response table
392 0 : pA1 = dAngleRad + respImageRotOffset(0).radian();
393 0 : pA2 = dAngleRad + respImageRotOffset(nAntTypes-1).radian();
394 :
395 0 : switch(polToDoIndex(iPol)){
396 0 : case 0: // XX
397 0 : fact1Index = fact2Index = 0;
398 0 : break;
399 0 : case 1: //XY
400 0 : fact1Index = 0;
401 0 : fact2Index = 3;
402 0 : break;
403 0 : case 2: //YX
404 0 : fact1Index = 3;
405 0 : fact2Index = 0;
406 0 : break;
407 0 : case 3: //YY
408 0 : fact1Index = fact2Index = 3;
409 0 : break;
410 : }
411 :
412 0 : if(pA1 != pA2){ // rotate individual factors before multiplication
413 :
414 : // rotate factor 1
415 0 : SynthesisUtils::rotateComplexArray(os, respByPol(0)(fact1Index), dCoord, fact1,
416 : pA1, "LINEAR",
417 : false); // don't modify dCoord
418 : // if necessary rotate factor 2
419 0 : if((nAntTypes-1)==0 && fact2Index==fact1Index){ // also implies that pA1==PA2
420 0 : fact2.assign(fact1);
421 : }
422 : else{
423 0 : SynthesisUtils::rotateComplexArray(os, respByPol(nAntTypes-1)(fact2Index), dCoord, fact2,
424 : pA2, "LINEAR",
425 : false); // don't modify dCoord
426 : }
427 : }
428 : else{ // rotate PB later
429 0 : fact1.assign(respByPol(0)(fact1Index));
430 0 : fact2.assign(respByPol(nAntTypes-1)(fact2Index));
431 : }
432 :
433 : // multiply EFPs (equivalent to convolution of AIFs) to get primary beam
434 0 : if(doSquint){
435 0 : pB = fact1 * conj(fact2);
436 : }
437 : else{
438 0 : pB = abs(fact1) * abs(fact2);
439 : }
440 :
441 : // now have the primary beam for polarization iPol in pB
442 :
443 : // combine all PBs into one image
444 0 : IPosition pos(rNDimFinal,0);
445 0 : pos(rAxis) = iPol;
446 :
447 0 : if(pA1 == pA2){ // still need to rotate pB by PA
448 :
449 0 : Array<Complex> pBrot( pB.shape() );
450 0 : SynthesisUtils::rotateComplexArray(os, pB, dCoord, pBrot,
451 : pA1, "LINEAR",
452 : false); // don't modify dCoord
453 0 : nearFinal.putSlice(pBrot, pos);
454 :
455 0 : }
456 : else{ // pB was already rotated above
457 0 : nearFinal.putSlice(pB, pos);
458 : }
459 :
460 0 : }
461 :
462 : // Then regrid it to the image
463 :
464 : // The following mess is necessary since ImageRegrid does not work for Complex images
465 : // (Interpolate2D is not templated)
466 :
467 0 : Array<Complex> nearFinalArray = nearFinal.get();
468 :
469 0 : CoordinateSystem outCS(outImage.coordinates());
470 0 : Vector<Int> pixAxes=outCS.pixelAxes(outCS.findCoordinate(Coordinate::DIRECTION));
471 0 : IPosition whichOutPixelAxes(pixAxes);
472 :
473 : // get the world coordinates of the center of outImage
474 0 : Vector<Double> wCenterOut(2);
475 0 : Vector<Double> pCenterOut(2);
476 0 : pCenterOut(0) = outImage.shape()(whichOutPixelAxes(0))/2.-0.5;
477 0 : pCenterOut(1) = outImage.shape()(whichOutPixelAxes(1))/2.-0.5;
478 0 : Vector<String> wAU = outCS.directionCoordinate(0).worldAxisUnits();
479 0 : outCS.directionCoordinate(0).toWorld(wCenterOut, pCenterOut);
480 : //cout << "pixel center " << pCenterOut << " world center " << wCenterOut << " " << wAU(0) << endl;
481 :
482 0 : uInt dirCoordIndex = dCoordFinal.findCoordinate(Coordinate::DIRECTION);
483 :
484 : // convert direction coordinate to J2000 if necessary
485 0 : if(dCoordFinal.directionCoordinate(dirCoordIndex).directionType() != MDirection::J2000){
486 0 : Vector<Double> incrV = dCoordFinal.directionCoordinate(dirCoordIndex).increment();
487 0 : Vector<Double> refPV = dCoordFinal.directionCoordinate(dirCoordIndex).referencePixel();
488 0 : Vector<String> wAxUnitsV = dCoordFinal.directionCoordinate(dirCoordIndex).worldAxisUnits();
489 0 : Projection pproj = dCoordFinal.directionCoordinate(dirCoordIndex).projection();
490 0 : Matrix<Double> xxform = dCoordFinal.directionCoordinate(dirCoordIndex).linearTransform();
491 0 : Quantity incr0(incrV(0),wAxUnitsV(0));
492 0 : Quantity incr1(incrV(1),wAxUnitsV(1));
493 :
494 : DirectionCoordinate newDC(MDirection::J2000, pproj, 0., 0.,
495 0 : incr0.getValue(Unit("rad")), incr1.getValue(Unit("rad")),
496 0 : xxform, refPV(0), refPV(1), 0., 90.);
497 0 : if(!dCoordFinal.replaceCoordinate(newDC, dirCoordIndex)){
498 0 : ostringstream oss;
499 0 : oss << "Error: response image direction coordinate cannot be set to direction type J2000.";
500 0 : throw(SynthesisError(oss.str()));
501 0 : }
502 0 : }
503 :
504 0 : DirectionCoordinate dCoordFinalDir = dCoordFinal.directionCoordinate(0);
505 :
506 : // Scale response image with frequency
507 0 : Unit uHz("Hz");
508 0 : Double refFreqHz = MFrequency::Convert(refFreq, MFrequency::TOPO)().get(uHz).getValue();
509 0 : Double nomFreqHz = MFrequency::Convert(respImageNomFreq(0), MFrequency::TOPO)().get(uHz).getValue();
510 0 : if(refFreqHz<=0.){
511 0 : throw(SynthesisError("Reference frequency in input image is <=0."));
512 : }
513 0 : Double fScale = nomFreqHz/refFreqHz;
514 0 : Vector<Double> newIncr = dCoordFinalDir.increment();
515 0 : cout << "Image increment before freq scaling " << newIncr << ", after scaling " << newIncr * fScale << endl;
516 0 : dCoordFinalDir.setIncrement(newIncr * fScale);
517 :
518 : // and set the reference Value of inImage to the world center of the outImage
519 0 : Vector<String> wAUI = dCoordFinalDir.worldAxisUnits();
520 0 : if(!(wAU(0)==wAUI(0))){
521 0 : Unit uIn(wAU(0));
522 0 : Unit uOut(wAUI(0));
523 0 : Quantity q0(wCenterOut(0), uIn);
524 0 : Quantity q1(wCenterOut(1), uIn);
525 0 : wCenterOut(0) = q0.getValue(uOut);
526 0 : wCenterOut(1) = q1.getValue(uOut);
527 0 : }
528 :
529 0 : cout << "Image pixel center " << pCenterOut << ", world center " << wCenterOut << " " << wAUI(0) << endl;
530 0 : dCoordFinalDir.setReferenceValue(wCenterOut);
531 :
532 0 : dCoordFinal.replaceCoordinate(dCoordFinalDir, dirCoordIndex);
533 :
534 : //PagedImage<Float> im1(nearFinalArray.shape(),dCoordFinal, "inImageReal.im");
535 : //im1.copyData(LatticeExpr<Float>(real(ArrayLattice<Complex>(nearFinalArray))));
536 : //PagedImage<Float> im2(nearFinalArray.shape(),dCoordFinal, "inImageImag.im");
537 : //im2.copyData(LatticeExpr<Float>(imag(ArrayLattice<Complex>(nearFinalArray))));
538 :
539 0 : TempImage<Float> inImage(nearFinalArray.shape(),dCoordFinal);
540 :
541 0 : TempImage<Float> tOutImage(outImage.shape(), outCS);
542 0 : Array<Complex> outArray(outImage.shape(), Complex(0.,0.));
543 :
544 0 : ImageRegrid<Float> iR;
545 : //iR.showDebugInfo(1);
546 :
547 0 : inImage.copyData(LatticeExpr<Float>(real(ArrayLattice<Complex>(nearFinalArray))));
548 0 : tOutImage.set(0.0);
549 :
550 0 : iR.regrid(tOutImage, Interpolate2D::LINEAR, whichOutPixelAxes, inImage);
551 0 : setReal(outArray,tOutImage.get());
552 :
553 0 : inImage.copyData(LatticeExpr<Float>(imag(ArrayLattice<Complex>(nearFinalArray))));
554 0 : tOutImage.set(0.0);
555 :
556 0 : iR.regrid(tOutImage, Interpolate2D::LINEAR, whichOutPixelAxes, inImage);
557 0 : setImag(outArray,tOutImage.get());
558 :
559 0 : outImage.put(outArray);
560 :
561 : // tidy up
562 0 : for(uInt i=0; i<nAntTypes; i++){
563 0 : delete respImage_p(i);
564 : }
565 0 : respImage_p.resize(0);
566 :
567 0 : }
568 0 : }
569 :
570 0 : void ALMAAperture::applySky(ImageInterface<Float>& outImage,
571 : const VisBuffer& vb,
572 : const Bool doSquint,
573 : const Int& cfKey,
574 : const Bool raytrace)
575 : {
576 0 : if(getVisParams(vb)==-1 || raytrace){ // need to use ray tracing
577 0 : ALMACalcIlluminationConvFunc almaPB;
578 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
579 0 : almaPB.setMaximumCacheSize(cachesize);
580 0 : almaPB.applyPB(outImage, vb, doSquint, cfKey);
581 0 : }
582 : else{ // use canned antenna responses
583 0 : TempImage<Complex> tI(outImage.shape(), outImage.coordinates());
584 0 : tI.set(Complex(0.,0.));
585 0 : applySky(tI, vb, doSquint, cfKey, false);
586 0 : outImage.put(real(abs(tI.get())));
587 0 : }
588 0 : }
589 :
590 0 : Int ALMAAperture::makePBPolnCoords(const VisBuffer&vb,
591 : const Int& convSize,
592 : const Int& convSampling,
593 : const CoordinateSystem& skyCoord,
594 : const Int& skyNx, const Int& /*skyNy*/,
595 : CoordinateSystem& feedCoord)
596 : {
597 0 : feedCoord = skyCoord;
598 : //
599 : // Make a two dimensional image to calculate auto-correlation of
600 : // the ideal illumination pattern. We want this on a fine grid in
601 : // the UV plane
602 : //
603 0 : Int directionIndex=skyCoord.findCoordinate(Coordinate::DIRECTION);
604 0 : AlwaysAssert(directionIndex>=0, AipsError);
605 0 : DirectionCoordinate dc=skyCoord.directionCoordinate(directionIndex);
606 0 : Vector<Double> sampling;
607 0 : sampling = dc.increment();
608 0 : sampling*=Double(convSampling);
609 0 : sampling*=Double(skyNx)/Double(convSize);
610 0 : dc.setIncrement(sampling);
611 :
612 :
613 0 : Vector<Double> unitVec(2);
614 0 : unitVec=convSize/2;
615 0 : dc.setReferencePixel(unitVec);
616 :
617 : // Set the reference value to that of the image
618 0 : feedCoord.replaceCoordinate(dc, directionIndex);
619 :
620 : //
621 : // Make an image with linear polarization axis.
622 : //
623 0 : Int NPol=0,M,N=0;
624 0 : M=polMap_p.nelements();
625 0 : for(Int i=0;i<M;i++){
626 0 : if (polMap_p(i) > -1) NPol++;
627 : }
628 0 : Vector<Int> poln(NPol);
629 :
630 : Int index;
631 0 : Vector<Int> inStokes;
632 0 : index = feedCoord.findCoordinate(Coordinate::STOKES);
633 0 : inStokes = feedCoord.stokesCoordinate(index).stokes();
634 0 : N = 0;
635 : try{
636 0 : for(Int i=0;i<M;i++) if (polMap_p(i) > -1) {poln(N) = vb.corrType()(i);N++;}
637 0 : StokesCoordinate polnCoord(poln);
638 0 : Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
639 0 : feedCoord.replaceCoordinate(polnCoord,StokesIndex);
640 0 : }
641 0 : catch(AipsError& x){
642 0 : throw(SynthesisFTMachineError("Likely cause: Discrepancy between the poln. "
643 0 : "axis of the data and the image specifications."));
644 0 : }
645 :
646 0 : return NPol;
647 0 : }
648 :
649 :
650 0 : Vector<Int> ALMAAperture::vbRow2CFKeyMap(const VisBuffer& vb, Int& nUnique)
651 : {
652 : // return index to outputImages for each row in vb
653 0 : Vector<Int> map;
654 0 : map.resize(vb.nRow());
655 :
656 : //cout << "vb rows: " << vb.nRow() << endl;
657 :
658 0 : if(haveCannedResponses_p){
659 : // distinguish different antenna types
660 0 : if(antTypeMap_p.nelements()!=(uInt)vb.numberAnt()){
661 0 : antTypeMap_p.assign(antTypeMap(vb));
662 0 : cout << "initialising antTypeMap to " << antTypeMap_p << endl;
663 : }
664 :
665 0 : std::map<Int, Int > cFKeysEncountered;
666 0 : Int cfKeyCount = 0;
667 0 : for(uInt i=0; i<(uInt)vb.nRow(); i++){
668 0 : Int cfKey = cFKeyFromAntennaTypes(antTypeMap_p(vb.antenna1()(i)),
669 0 : antTypeMap_p(vb.antenna2()(i)));
670 0 : map(i) = cfKey;
671 0 : if(cFKeysEncountered.find(cfKey) == cFKeysEncountered.end( )){ // new cFKey
672 0 : cFKeysEncountered.insert(std::pair<Int, Int >(cfKey, cfKeyCount));
673 0 : cfKeyCount++;
674 : }
675 : }
676 0 : nUnique = cfKeyCount;
677 0 : }
678 : else{ // raytracing knows only one antenna type for the moment
679 0 : map=0;
680 0 : nUnique=1;
681 : }
682 0 : return map;
683 0 : }
684 :
685 0 : Vector<ALMAAntennaType> ALMAAperture::antTypeMap(const VisBuffer& vb){
686 0 : Vector<ALMAAntennaType> map(vb.numberAnt(),ALMA_INVALID);
687 0 : for(uInt i=0; i<map.nelements(); i++){
688 0 : map(i) = antTypeFromName(vb.msColumns().antenna().name()(i));
689 0 : cout << vb.msColumns().antenna().name()(i) << " " << map(i) << endl;
690 0 : if(map(i)==ALMA_INVALID){
691 0 : logIO() << LogOrigin("ALMAAperture", "antTypeMap")
692 : << LogIO::WARN
693 : << "Unrecognised antenna type for antenna \""
694 0 : << vb.msColumns().antenna().name()(i) << "\""
695 0 : << LogIO::POST;
696 : }
697 : }
698 0 : return map;
699 0 : }
700 :
701 0 : ALMAAntennaType ALMAAperture::antTypeFromName(const String& name){
702 0 : ALMAAntennaType rval = ALMA_INVALID;
703 0 : String typeN = name.substr(0,2);
704 0 : if(typeN=="DV") rval=ALMA_DV;
705 0 : else if(typeN=="DA") rval=ALMA_DA;
706 0 : else if(typeN=="PM") rval=ALMA_PM;
707 0 : else if(typeN=="CM") rval=ALMA_CM;
708 0 : return rval;
709 0 : }
710 :
711 0 : String ALMAAperture::antTypeStrFromType(const ALMAAntennaType& aType){
712 0 : String tS;
713 0 : switch(aType){
714 0 : case ALMA_DV:
715 0 : tS = "DV";
716 0 : break;
717 0 : case ALMA_DA:
718 0 : tS = "DA";
719 0 : break;
720 0 : case ALMA_PM:
721 0 : tS = "PM";
722 0 : break;
723 0 : case ALMA_CM:
724 0 : tS = "CM";
725 0 : break;
726 0 : case ALMA_INVALID:
727 : default:
728 0 : tS = "UNKOWN";
729 0 : break;
730 : }
731 0 : return tS;
732 0 : }
733 :
734 :
735 0 : Int ALMAAperture::cFKeyFromAntennaTypes(const ALMAAntennaType aT1, const ALMAAntennaType aT2){
736 0 : if(orderMattersInCFKey){
737 0 : return (Int)aT2+1 + 10000*((Int)aT1+1);
738 : }
739 : else{
740 0 : return min((Int)aT1+1, (Int)aT2+1) + 10000*max((Int)aT1+1, (Int)aT2+1); // assume order doesn't matter
741 : }
742 : }
743 :
744 0 : Vector<ALMAAntennaType> ALMAAperture::antennaTypesFromCFKey(const Int& cFKey){
745 0 : Int t1 = (Int) floor(cFKey/10000.);
746 0 : Int t2 = cFKey - 10000*t1 - 1;
747 0 : t1--;
748 0 : Vector<ALMAAntennaType> rval(1);
749 0 : if(t1==t2){
750 0 : rval(0) = (ALMAAntennaType) t1;
751 : }
752 : else{
753 0 : rval.resize(2);
754 0 : if(orderMattersInCFKey){
755 0 : rval(0) = (ALMAAntennaType)t1;
756 0 : rval(1) = (ALMAAntennaType)t2;
757 : }
758 : else{
759 0 : rval(0) = (ALMAAntennaType) min(t1,t2);
760 0 : rval(1) = (ALMAAntennaType) max(t1,t2);
761 : }
762 : }
763 0 : return rval;
764 0 : }
765 :
766 0 : Vector<ALMAAntennaType> ALMAAperture::antTypeList(const VisBuffer& vb){
767 0 : Vector<ALMAAntennaType> aTypeMap = antTypeMap(vb);
768 0 : Vector<Bool> encountered(ALMA_numAntTypes, false);
769 0 : uInt typeCount = 0;
770 0 : for(uInt i=0; i<aTypeMap.nelements(); i++){
771 0 : uInt index = (uInt)aTypeMap(i);
772 0 : if(!encountered(index)){
773 0 : encountered(index) = true;
774 0 : typeCount++;
775 : }
776 : }
777 0 : Vector<ALMAAntennaType> aTypeList(typeCount);
778 0 : uInt index = 0;
779 0 : for(uInt i=0; i<ALMA_numAntTypes; i++){
780 0 : if(encountered(i)){
781 0 : aTypeList(index) = (ALMAAntennaType)i;
782 0 : index++;
783 : }
784 : }
785 0 : return aTypeList;
786 0 : }
787 :
788 :
789 : };
|
