Line data Source code
1 : // -*- C++ -*-
2 : //# PBMosaicFT.cc: Implementation of PBMosaicFT 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/VisibilityIterator.h>
30 : #include <casacore/casa/Quanta/UnitMap.h>
31 : #include <casacore/casa/Quanta/MVTime.h>
32 : #include <casacore/casa/Quanta/UnitVal.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/Projection.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/ms/MeasurementSets/MSRange.h>
41 : #include <casacore/ms/MSSel/MSSelection.h>
42 : #include <casacore/casa/BasicSL/Constants.h>
43 : #include <casacore/scimath/Mathematics/FFTServer.h>
44 : #include <synthesis/MeasurementComponents/PBMosaicFT.h>
45 : #include <casacore/scimath/Mathematics/RigidVector.h>
46 : #include <msvis/MSVis/StokesVector.h>
47 : #include <synthesis/TransformMachines/StokesImageUtil.h>
48 : #include <msvis/MSVis/VisBuffer.h>
49 : #include <msvis/MSVis/VisSet.h>
50 : #include <casacore/images/Images/ImageInterface.h>
51 : #include <casacore/images/Images/PagedImage.h>
52 : #include <casacore/casa/Containers/Block.h>
53 : #include <casacore/casa/Containers/Record.h>
54 : #include <casacore/casa/Arrays/ArrayLogical.h>
55 : #include <casacore/casa/Arrays/ArrayMath.h>
56 : #include <casacore/casa/Arrays/Array.h>
57 : #include <casacore/casa/Arrays/MaskedArray.h>
58 : #include <casacore/casa/Arrays/Vector.h>
59 : #include <casacore/casa/Arrays/Slice.h>
60 : #include <casacore/casa/Arrays/Matrix.h>
61 : #include <casacore/casa/Arrays/Cube.h>
62 : #include <casacore/casa/Arrays/MatrixIter.h>
63 : #include <casacore/casa/BasicSL/String.h>
64 : #include <casacore/casa/Utilities/Assert.h>
65 : #include <casacore/casa/Exceptions/Error.h>
66 : #include <casacore/lattices/Lattices/ArrayLattice.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LRegions/LCBox.h>
69 : #include <casacore/lattices/LEL/LatticeExpr.h>
70 : #include <casacore/lattices/Lattices/LatticeCache.h>
71 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
72 : #include <casacore/lattices/Lattices/LatticeIterator.h>
73 : #include <casacore/lattices/Lattices/LatticeStepper.h>
74 : #include <casacore/casa/Utilities/CompositeNumber.h>
75 : #include <casacore/casa/OS/Timer.h>
76 : #include <casacore/casa/OS/HostInfo.h>
77 : #include <sstream>
78 : #include <casacore/images/Regions/WCBox.h>
79 : #include <casacore/images/Images/SubImage.h>
80 : #include <casacore/images/Regions/ImageRegion.h>
81 : #include <casacore/images/Images/ImageSummary.h>
82 : #include <casacore/casa/BasicMath/Random.h>
83 :
84 : #include <synthesis/MeasurementComponents/PBMosaicFT.h>
85 : //#include <synthesis/MeasurementComponents/GlobalFTMachineCallbacks.h>
86 : #include <casacore/casa/System/ProgressMeter.h>
87 :
88 : #define DELTAPA 1.0
89 : #define MAGICPAVALUE -999.0
90 : #define CONVSIZE (1024*4)
91 : #define OVERSAMPLING 20
92 : #define USETABLES 0 // If equal to 1, use tabulated exp() and
93 : // complex exp() functions.
94 : #define MAXPOINTINGERROR 250.0 // Max. pointing error in arcsec used to
95 : // determine the resolution of the
96 : // tabulated exp() function.
97 : using namespace casacore;
98 : namespace casa { //# NAMESPACE CASA - BEGIN
99 : //
100 : //---------------------------------------------------------------
101 : //
102 : #define FUNC(a) (sqrt((a)))
103 0 : PBMosaicFT::PBMosaicFT(MeasurementSet& /*ms*/,
104 : Int nWPlanes, Long icachesize,
105 : String& cfCacheDirName,
106 : Bool applyPointingOffset,
107 : Bool doPBCorr,
108 : Int itilesize,
109 : Float paSteps,
110 : Float pbLimit,
111 0 : Bool usezero)
112 : ////: nPBWProjectFT(ms,nWPlanes,icachesize,cfCacheDirName,applyPointingOffset,doPBCorr,itilesize,paSteps,pbLimit,usezero),
113 : : nPBWProjectFT(nWPlanes,icachesize,cfCacheDirName,applyPointingOffset,doPBCorr,itilesize,paSteps,pbLimit,usezero),
114 0 : fieldIds_p(0)
115 : {
116 : //
117 : // Set the function pointer for FORTRAN call for GCF services.
118 : // This is a pure C function pointer which FORTRAN can call.
119 : // The function itself uses GCF object services.
120 : //
121 0 : epJ=NULL; // We do not yet support antenna pointing error
122 : // handling in this FTMachine.
123 0 : convSize=0;
124 0 : tangentSpecified_p=false;
125 0 : lastIndex_p=0;
126 0 : paChangeDetector.reset();
127 0 : pbLimit_p=pbLimit;
128 : //
129 : // Get various parameters from the visibilities.
130 : //
131 : //bandID_p = getVisParams();
132 0 : if (applyPointingOffset) doPointing=1; else doPointing=0;
133 :
134 0 : convFuncCacheReady=false;
135 0 : PAIndex = -1;
136 0 : maxConvSupport=-1;
137 : //
138 : // Set up the Conv. Func. disk cache manager object.
139 : //
140 0 : cfCache.setCacheDir(cfCacheDirName.data());
141 0 : cfCache.initCache();
142 0 : convSampling=OVERSAMPLING;
143 0 : convSize=CONVSIZE;
144 : //use memory size defined in aipsrc if exists
145 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
146 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
147 0 : if (cachesize > hostRAM) cachesize=hostRAM;
148 :
149 : // MSFieldColumns msfc(ms.field());
150 : // MSRange msr(ms);
151 : // //
152 : // // An array of shape [2,1,1]!
153 : // //
154 : // fieldIds_p = msr.range(MSS::FIELD_ID).asArrayInt(RecordFieldId(0));
155 0 : nApertures = 0;
156 0 : lastPAUsedForWtImg = MAGICPAVALUE;
157 0 : }
158 : //
159 : //---------------------------------------------------------------
160 : //
161 0 : PBMosaicFT::PBMosaicFT(const RecordInterface& stateRec)
162 0 : : nPBWProjectFT(stateRec)
163 : {
164 : // Construct from the input state record
165 0 : String error;
166 :
167 0 : if (!fromRecord(error, stateRec)) {
168 0 : throw (AipsError("Failed to create PBMosaicFT: " + error));
169 : };
170 : // bandID_p = getVisParams();
171 0 : PAIndex = -1;
172 0 : maxConvSupport=-1;
173 0 : convSampling=OVERSAMPLING;
174 0 : convSize=CONVSIZE;
175 0 : nApertures = 0;
176 0 : }
177 : //
178 : //---------------------------------------------------------------
179 : //
180 0 : PBMosaicFT& PBMosaicFT::operator=(const PBMosaicFT& other)
181 : {
182 0 : if(this!=&other) {
183 0 : nPBWProjectFT::operator=(other);
184 0 : padding_p=other.padding_p;
185 : // ms_p=other.ms_p;
186 0 : nWPlanes_p=other.nWPlanes_p;
187 0 : imageCache=other.imageCache;
188 0 : cachesize=other.cachesize;
189 0 : tilesize=other.tilesize;
190 0 : gridder=other.gridder;
191 0 : isTiled=other.isTiled;
192 0 : lattice=other.lattice;
193 0 : arrayLattice=other.arrayLattice;
194 0 : maxAbsData=other.maxAbsData;
195 0 : centerLoc=other.centerLoc;
196 0 : offsetLoc=other.offsetLoc;
197 0 : mspc=other.mspc;
198 0 : msac=other.msac;
199 0 : pointingToImage=other.pointingToImage;
200 0 : usezero_p=other.usezero_p;
201 0 : doPBCorrection = other.doPBCorrection;
202 0 : maxConvSupport= other.maxConvSupport;
203 0 : nApertures = other.nApertures;
204 : };
205 0 : return *this;
206 : };
207 : //
208 : //----------------------------------------------------------------------
209 : //
210 : // PBMosaicFT::PBMosaicFT(const PBMosaicFT& other)
211 : // {
212 : // operator=(other);
213 : // }
214 : //
215 : //---------------------------------------------------------------
216 : //
217 0 : Int PBMosaicFT::findPointingOffsets(const VisBuffer& vb,
218 : Array<Float> &l_off,
219 : Array<Float> &m_off,
220 : Bool Evaluate)
221 : {
222 : //
223 : // Get the antenna pointing errors, if any.
224 0 : Int NAnt=0;
225 : //
226 : // This will return 0 if EPJ Table is not given. Otherwise will
227 : // return the number of antennas it detected (from the EPJ table)
228 : // and the offsets in l_off and m_off.
229 : //
230 0 : NAnt = nPBWProjectFT::findPointingOffsets(vb,l_off,m_off,Evaluate);
231 : // NAnt = l_off.shape()(2);
232 :
233 : // Resize the offset arrays if no pointing table was given.
234 : //
235 0 : if (NAnt <=0 )
236 : {
237 0 : NAnt=vb.msColumns().antenna().nrow();
238 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
239 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
240 0 : l_off = m_off = 0.0;
241 : }
242 : //
243 : // Add field offsets to the pointing errors.
244 : //
245 :
246 : // Float dayHr=2*3.141592653589793116;
247 : // MVDirection ref(directionCoord.referenceValue()(0),
248 : // directionCoord.referenceValue()(1)),
249 : // vbDir(vb.direction1()(0).getAngle().getValue()(0),
250 : // vb.direction1()(0).getAngle().getValue()(1));
251 :
252 : // if (0)
253 : // {
254 : // l_off = l_off - (Float)(directionCoord.referenceValue()(0) -
255 : // dayHr-vb.direction1()(0).getAngle().getValue()(0));
256 : // m_off = m_off - (Float)(directionCoord.referenceValue()(1) -
257 : // vb.direction1()(0).getAngle().getValue()(1));
258 : // }
259 :
260 : //
261 : // Convert the direction from image co-ordinate system and the VB
262 : // to MDirection. Then convert the MDirection to Quantity so that
263 : // arithematic operation (subtraction) can be done. Then use the
264 : // subtracted Quantity to construct another MDirection and use
265 : // *it's* getAngle().getValue() to extract the difference in the
266 : // sky direction (from image co-ordinate system) and the phase
267 : // center direction of the VB in radians!
268 : //
269 : // If only VisBuffer, and DirectionCoordinate class could return
270 : // MDirection, and MDirection class had operator-(), the code
271 : // below could look more readable as:
272 : // MDirection diff=vb.mDirection()-directionCoord.mDirection();
273 : //
274 0 : MDirection vbMDir(vb.direction1()(0)),skyMDir, diff;
275 0 : directionCoord.toWorld(skyMDir, directionCoord.referencePixel());
276 0 : diff = MDirection(skyMDir.getAngle()-vbMDir.getAngle());
277 0 : l_off = l_off - (Float)diff.getAngle().getValue()(0);
278 0 : m_off = m_off - (Float)diff.getAngle().getValue()(1);
279 :
280 : static int firstPass=0;
281 : // static int fieldID=-1;
282 0 : static Vector<Float> offsets0,offsets1;
283 0 : if (firstPass==0)
284 : {
285 0 : offsets0.resize(NAnt);
286 0 : offsets1.resize(NAnt);
287 : // MLCG mlcg((Int)(vb.time()(0)));
288 : // Normal nrand(&mlcg,0.0,10.0);
289 : // for(Int i=0;i<NAnt;i++) offsets0(i) = (Float)(nrand());
290 : // for(Int i=0;i<NAnt;i++) offsets1(i) = (Float)(nrand());
291 0 : offsets0 = offsets1 = 0.0;
292 : }
293 0 : for(Int i=0;i<NAnt;i++)
294 : {
295 0 : l_off(IPosition(3,0,0,i)) = l_off(IPosition(3,0,0,i)) + offsets0(i)/2.062642e+05;
296 0 : m_off(IPosition(3,0,0,i)) = m_off(IPosition(3,0,0,i)) + offsets1(i)/2.062642e+05;
297 : }
298 :
299 : //m_off=l_off=0.0;
300 : // if (fieldID != vb.fieldId())
301 : // {
302 : // fieldID = vb.fieldId();
303 : // cout << l_off*2.062642e5 << endl;
304 : // cout << m_off*2.062642e5 << endl;
305 : // }
306 0 : if (firstPass==0)
307 : {
308 : // cout << (Float)(directionCoord.referenceValue()(0)) << " "
309 : // << (Float)(directionCoord.referenceValue()(1)) << endl;
310 : // cout << vb.direction1()(0).getAngle().getValue()(0) << " "
311 : // << vb.direction1()(0).getAngle().getValue()(1) << endl;
312 : // cout << l_off << endl;
313 : // cout << m_off << endl;
314 : }
315 0 : firstPass++;
316 0 : return NAnt;
317 : }
318 : //
319 : //---------------------------------------------------------------
320 : //
321 0 : void PBMosaicFT::normalizeAvgPB()
322 : {
323 : // We accumulated normalized PBs. So don't normalize the average
324 : // PB.
325 0 : pbNormalized = false;
326 :
327 0 : }
328 : //
329 : //---------------------------------------------------------------
330 : //
331 0 : Bool PBMosaicFT::makeAveragePB0(const VisBuffer& /*vb*/,
332 : const ImageInterface<Complex>& image,
333 : Int& /*polInUse*/,
334 : TempImage<Float>& theavgPB)
335 : {
336 0 : Bool pbMade=false;
337 0 : if (!resetPBs) return pbMade;
338 : //
339 : // If this is the first time, resize the average PB
340 : //
341 0 : if (resetPBs)
342 : {
343 0 : theavgPB.resize(image.shape());
344 0 : theavgPB.setCoordinateInfo(image.coordinates());
345 0 : theavgPB.set(0.0);
346 0 : noOfPASteps = 1;
347 0 : pbPeaks.resize(theavgPB.shape()(2));
348 0 : pbPeaks.set(0.0);
349 0 : resetPBs=false;
350 0 : pbNormalized=false;
351 : }
352 0 : return pbMade;
353 : }
354 : //
355 : //--------------------------------------------------------------------------------
356 : //
357 0 : void PBMosaicFT::normalizePB(ImageInterface<Float>& /*pb*/, const Float& /*peakValue*/)
358 : {
359 0 : }
360 : //
361 : //---------------------------------------------------------------
362 : //
363 : // Finalize the FFT to the Sky. Here we actually do the FFT and
364 : // return the resulting image
365 : //
366 : // A specialization exists here only to not normalize by the avgPB.
367 : //
368 0 : ImageInterface<Complex>& PBMosaicFT::getImage(Matrix<Float>& weights,
369 : Bool normalize)
370 : {
371 : //AlwaysAssert(lattice, AipsError);
372 0 : AlwaysAssert(image, AipsError);
373 :
374 0 : logIO() << "#########getimage########" << LogIO::DEBUGGING << LogIO::POST;
375 :
376 0 : logIO() << LogOrigin("PBMosaicFT", "getImage") << LogIO::NORMAL;
377 :
378 0 : weights.resize(sumWeight.shape());
379 :
380 0 : convertArray(weights, sumWeight);
381 : // cerr << "Sum Wt = " << sumWeight << endl;
382 : //
383 : // If the weights are all zero then we cannot normalize otherwise
384 : // we don't care.
385 : //
386 0 : if(max(weights)==0.0)
387 : {
388 0 : if(normalize) logIO() << LogIO::SEVERE
389 : << "No useful data in PBMosaicFT: weights all zero"
390 0 : << LogIO::POST;
391 : else logIO() << LogIO::WARN << "No useful data in PBMosaicFT: weights all zero"
392 0 : << LogIO::POST;
393 : }
394 : else
395 : {
396 0 : const IPosition latticeShape = lattice->shape();
397 :
398 : logIO() << LogIO::DEBUGGING
399 0 : << "Starting FFT and scaling of image" << LogIO::POST;
400 : //
401 : // x and y transforms (lattice has the gridded vis. Make the
402 : // dirty images)
403 : //
404 : // if (!makingPSF)
405 : // {
406 : // String name("griddedVis.im");
407 : // image->put(griddedData);
408 : // storeImg(name,*image);
409 : // }
410 0 : LatticeFFT::cfft2d(*lattice,false);
411 0 : if (!avgPBReady)
412 : {
413 0 : avgPBReady=true;
414 0 : avgPB.resize(griddedWeights.shape());
415 0 : avgPB.setCoordinateInfo(griddedWeights.coordinates());
416 : // pbNormalized=true;
417 : // {
418 : // String name("cpb.im");
419 : // storeImg(name,griddedWeights);
420 : // }
421 0 : makeSensitivityImage(griddedWeights, avgPB, weights, true);
422 : //
423 : // For effeciency reasons, weight functions are
424 : // accumulated only once per channel and polarization per
425 : // VB (i.e. for a given VB, if the weight functions are
426 : // accumulated for every channel and polarization, they
427 : // are not accumulated for the rest of VB). This makes
428 : // the sum of weights for wegith functions different from
429 : // sum of weights for the visibility data and the
430 : // normalzation in makeSensitivityImage() will be
431 : // incorrect. For now, normalize the peak of the average
432 : // weight function (wavgPB) to 1.0.
433 : //
434 0 : pbNormalized=false;
435 0 : nPBWProjectFT::normalizeAvgPB();
436 0 : resetPBs=false;
437 0 : cfCache.finalize(avgPB);
438 : }
439 :
440 : //
441 : // Apply the gridding correction
442 : //
443 : {
444 : // normalizeAvgPB();
445 0 : Int inx = lattice->shape()(0);
446 0 : Int iny = lattice->shape()(1);
447 0 : Vector<Complex> correction(inx);
448 :
449 0 : Vector<Float> sincConv(nx);
450 0 : Float centerX=nx/2+1;
451 0 : for (Int ix=0;ix<nx;ix++)
452 : {
453 0 : Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
454 0 : if(ix==centerX) sincConv(ix)=1.0;
455 0 : else sincConv(ix)=sin(x)/x;
456 : }
457 :
458 0 : IPosition cursorShape(4, inx, 1, 1, 1);
459 0 : IPosition axisPath(4, 0, 1, 2, 3);
460 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
461 0 : LatticeIterator<Complex> lix(*lattice, lsx);
462 :
463 0 : LatticeStepper lavgpb(avgPB.shape(),cursorShape,axisPath);
464 0 : LatticeIterator<Float> liavgpb(avgPB, lavgpb);
465 0 : Float peakAvgPB = max(avgPB.get());
466 0 : for(lix.reset(),liavgpb.reset();
467 0 : !lix.atEnd();
468 0 : lix++,liavgpb++)
469 : {
470 0 : Int pol=lix.position()(2);
471 0 : Int chan=lix.position()(3);
472 : // if(weights(pol, chan)>0.0)
473 : {
474 0 : Int iy=lix.position()(1);
475 0 : gridder->correctX1D(correction,iy);
476 :
477 0 : Vector<Float> PBCorrection(liavgpb.rwVectorCursor().shape()),
478 0 : avgPBVec(liavgpb.rwVectorCursor().shape());
479 :
480 0 : PBCorrection = liavgpb.rwVectorCursor();
481 0 : avgPBVec = liavgpb.rwVectorCursor();
482 :
483 0 : sincConv=1.0;
484 0 : if (!doPBCorrection) avgPBVec=1.0;
485 : // avgPBVec=1.0;
486 0 : for(int i=0;i<PBCorrection.shape();i++)
487 : {
488 : //
489 : // This with PS functions
490 : //
491 : // PBCorrection(i)=pbFunc(avgPBVec(i))*sincConv(i)*sincConv(iy);
492 : // // PBCorrection(i)=(avgPBVec(i))*sincConv(i)*sincConv(iy);
493 : // if ((abs(PBCorrection(i)*correction(i))) >= pbLimit_p)
494 : // lix.rwVectorCursor()(i) /= PBCorrection(i)*correction(i); //*pbNorm
495 : // else
496 : // lix.rwVectorCursor()(i) /= correction(i)*sincConv(i)*sincConv(iy);//pbNorm;
497 : //
498 : // This without the PS functions
499 : //
500 0 : PBCorrection(i)=pbFunc(avgPBVec(i))*sincConv(i)*sincConv(iy);
501 0 : if ((abs(PBCorrection(i))) >= pbLimit_p*peakAvgPB)
502 0 : lix.rwVectorCursor()(i) /= PBCorrection(i);
503 0 : else if (!makingPSF)
504 0 : lix.rwVectorCursor()(i) /= sincConv(i)*sincConv(iy);
505 : }
506 :
507 0 : if(normalize)
508 : {
509 0 : if(weights(pol, chan)>0.0)
510 : {
511 0 : Complex rnorm(Float(inx)*Float(iny)/(weights(pol,chan)));
512 0 : lix.rwCursor()*=rnorm;
513 : }
514 : else
515 0 : lix.woCursor()=0.0;
516 : }
517 : else
518 : {
519 0 : Complex rnorm(Float(inx)*Float(iny));
520 0 : lix.rwCursor()*=rnorm;
521 : }
522 : }
523 : }
524 : }
525 :
526 0 : if(!isTiled)
527 : {
528 : //
529 : // Check the section from the image BEFORE converting to a lattice
530 : //
531 0 : IPosition blc(4, (nx-image->shape()(0))/2,
532 0 : (ny-image->shape()(1))/2, 0, 0);
533 0 : IPosition stride(4, 1);
534 0 : IPosition trc(blc+image->shape()-stride);
535 : //
536 : // Do the copy
537 : //
538 0 : image->put(griddedData(blc, trc));
539 : //if(arrayLattice) delete arrayLattice; arrayLattice=0;
540 0 : griddedData.resize(IPosition(1,0));
541 : }
542 : }
543 : // if (!makingPSF)
544 : // {
545 : // String name("cdirty.im");
546 : // image->put(griddedData);
547 : // storeImg(name,*image);
548 : // }
549 :
550 :
551 0 : return *image;
552 : }
553 :
554 :
555 : #define NEED_UNDERSCORES
556 : #if defined(NEED_UNDERSCORES)
557 : #define gpbmos gpbmos_
558 : #define dpbmos dpbmos_
559 : #define dpbmosgrad dpbmosgrad_
560 : #define dpbwgrad dpbwgrad_
561 : #endif
562 :
563 : extern "C" {
564 : void gpbmos(Double *uvw,
565 : Double *dphase,
566 : const Complex *values,
567 : Int *nvispol,
568 : Int *nvischan,
569 : Int *dopsf,
570 : const Int *flag,
571 : const Int *rflag,
572 : const Float *weight,
573 : Int *nrow,
574 : Int *rownum,
575 : Double *scale,
576 : Double *offset,
577 : Complex *grid,
578 : Int *nx,
579 : Int *ny,
580 : Int *npol,
581 : Int *nchan,
582 : const Double *freq,
583 : const Double *c,
584 : Int *support,
585 : Int *convsize,
586 : Int *convwtsize,
587 : Int *sampling,
588 : Int *wconvsize,
589 : Complex *convfunc,
590 : Int *chanmap,
591 : Int *polmap,
592 : Int *polused,
593 : Double *sumwt,
594 : Int *ant1,
595 : Int *ant2,
596 : Int *nant,
597 : Int *scanno,
598 : Double *sigma,
599 : Float *raoff,
600 : Float *decoff,
601 : Double *area,
602 : Int *doGrad,
603 : Int *doPointingCorrection,
604 : Int *nPA,
605 : Int *paIndex,
606 : Int *CFMap,
607 : Int *ConjCFMap,
608 : Double *currentCFPA, Double *actualPA,
609 : Int *avgPBReady,
610 : Complex *avgPB, Double *cfRefFreq_p,
611 : Complex *convWeights,
612 : Int *convWtSupport);
613 : void dpbmos(Double *uvw,
614 : Double *dphase,
615 : Complex *values,
616 : Int *nvispol,
617 : Int *nvischan,
618 : const Int *flag,
619 : const Int *rflag,
620 : Int *nrow,
621 : Int *rownum,
622 : Double *scale,
623 : Double *offset,
624 : const Complex *grid,
625 : Int *nx,
626 : Int *ny,
627 : Int *npol,
628 : Int *nchan,
629 : const Double *freq,
630 : const Double *c,
631 : Int *support,
632 : Int *convsize,
633 : Int *sampling,
634 : Int *wconvsize,
635 : Complex *convfunc,
636 : Int *chanmap,
637 : Int *polmap,
638 : Int *polused,
639 : Int *ant1,
640 : Int *ant2,
641 : Int *nant,
642 : Int *scanno,
643 : Double *sigma,
644 : Float *raoff, Float *decoff,
645 : Double *area,
646 : Int *dograd,
647 : Int *doPointingCorrection,
648 : Int *nPA,
649 : Int *paIndex,
650 : Int *CFMap,
651 : Int *ConjCFMap,
652 : Double *currentCFPA, Double *actualPA, Double *cfRefFreq_p);
653 : void dpbmosgrad(Double *uvw,
654 : Double *dphase,
655 : Complex *values,
656 : Int *nvispol,
657 : Int *nvischan,
658 : Complex *gazvalues,
659 : Complex *gelvalues,
660 : Int *doconj,
661 : const Int *flag,
662 : const Int *rflag,
663 : Int *nrow,
664 : Int *rownum,
665 : Double *scale,
666 : Double *offset,
667 : const Complex *grid,
668 : Int *nx,
669 : Int *ny,
670 : Int *npol,
671 : Int *nchan,
672 : const Double *freq,
673 : const Double *c,
674 : Int *support,
675 : Int *convsize,
676 : Int *sampling,
677 : Int *wconvsize,
678 : Complex *convfunc,
679 : Int *chanmap,
680 : Int *polmap,
681 : Int *polused,
682 : Int *ant1,
683 : Int *ant2,
684 : Int *nant,
685 : Int *scanno,
686 : Double *sigma,
687 : Float *raoff, Float *decoff,
688 : Double *area,
689 : Int *dograd,
690 : Int *doPointingCorrection,
691 : Int *nPA,
692 : Int *paIndex,
693 : Int *CFMap,
694 : Int *ConjCFMap,
695 : Double *currentCFPA, Double *actualPA,Double *cfRefFreq_p);
696 : }
697 : //
698 : //----------------------------------------------------------------------
699 : //
700 0 : void PBMosaicFT::runFortranGet(Matrix<Double>& uvw,Vector<Double>& dphase,
701 : Cube<Complex>& visdata,
702 : IPosition& s,
703 : Int& /*Conj*/,
704 : Cube<Int>& flags,Vector<Int>& rowFlags,
705 : Int& rownr,Vector<Double>& actualOffset,
706 : Array<Complex>* dataPtr,
707 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
708 : VisBuffer& vb,Int& Nant_p, Int& scanNo,
709 : Double& sigma,
710 : Array<Float>& l_off,
711 : Array<Float>& m_off,
712 : Double area,
713 : Int& doGrad,
714 : Int paIndex)
715 : {
716 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
717 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
718 : CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
719 0 : Vector<Bool> deleteThem(21);
720 :
721 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
722 : Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
723 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
724 : *ConjCFMap_p, *CFMap_p;
725 : Float *l_off_p, *m_off_p;
726 : Double actualPA;
727 :
728 0 : Vector<Int> ConjCFMap, CFMap;
729 : Int N;
730 0 : actualPA = getVBPA(vb);
731 :
732 0 : N=polMap.nelements();
733 0 : CFMap = polMap; ConjCFMap = polMap;
734 0 : for(Int i=0;i<N;i++) CFMap[i] = polMap[N-i-1];
735 :
736 0 : Array<Complex> rotatedConvFunc;
737 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
738 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
739 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
740 : rotatedConvFunc,0.0);
741 :
742 0 : ConjCFMap = polMap;
743 0 : makeCFPolMap(vb,CFMap);
744 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
745 :
746 :
747 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
748 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
749 :
750 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
751 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
752 0 : visdata_p = visdata.getStorage(deleteThem(VISDATA));
753 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
754 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
755 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
756 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
757 0 : dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
758 0 : vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
759 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
760 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
761 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
762 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
763 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
764 0 : vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
765 0 : vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
766 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
767 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
768 :
769 : // static int ttt=0;
770 : // if (ttt==0) cout << l_off(IPosition(3,0,0,0)) << " " << m_off(IPosition(3,0,0,0)) << endl;
771 : // ttt++;
772 :
773 0 : Int npa=convSupport.shape()(2),actualConvSize;
774 0 : Int paIndex_Fortran = paIndex;
775 0 : actualConvSize = convFunc.shape()(0);
776 :
777 0 : IPosition shp=convSupport.shape();
778 0 : Int alwaysDoPointing=1;
779 0 : alwaysDoPointing=doPointing;
780 0 : dpbmos(uvw_p,
781 : dphase_p,
782 : visdata_p,
783 0 : &s.asVector()(0),
784 0 : &s.asVector()(1),
785 : flags_p,
786 : rowFlags_p,
787 0 : &s.asVector()(2),
788 : &rownr,
789 : uvScale_p,
790 : actualOffset_p,
791 : dataPtr_p,
792 : &aNx,
793 : &aNy,
794 : &npol,
795 : &nchan,
796 : vb_freq_p,
797 : &C::c,
798 : convSupport_p,
799 : &actualConvSize,
800 : &convSampling,
801 : &wConvSize,
802 : f_convFunc_p,
803 : chanMap_p,
804 : polMap_p,
805 : &polInUse,
806 : vb_ant1_p,
807 : vb_ant2_p,
808 : &Nant_p,
809 : &scanNo,
810 : &sigma,
811 : l_off_p, m_off_p,
812 : &area,
813 : &doGrad,
814 : &alwaysDoPointing,
815 : &npa,
816 : &paIndex_Fortran,
817 : CFMap_p,
818 : ConjCFMap_p,
819 : ¤tCFPA
820 : ,&actualPA, &cfRefFreq_p
821 : );
822 :
823 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
824 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
825 :
826 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
827 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
828 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
829 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
830 0 : visdata.putStorage(visdata_p,deleteThem(VISDATA));
831 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
832 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
833 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
834 0 : dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
835 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
836 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
837 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
838 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
839 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
840 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
841 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
842 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
843 0 : }
844 : //
845 : //----------------------------------------------------------------------
846 : //
847 0 : void PBMosaicFT::runFortranGetGrad(Matrix<Double>& uvw,Vector<Double>& dphase,
848 : Cube<Complex>& visdata,
849 : IPosition& s,
850 : Cube<Complex>& gradVisAzData,
851 : Cube<Complex>& gradVisElData,
852 : Int& Conj,
853 : Cube<Int>& flags,Vector<Int>& rowFlags,
854 : Int& rownr,Vector<Double>& actualOffset,
855 : Array<Complex>* dataPtr,
856 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
857 : VisBuffer& vb,Int& Nant_p, Int& scanNo,
858 : Double& sigma,
859 : Array<Float>& l_off,
860 : Array<Float>& m_off,
861 : Double area,
862 : Int& doGrad,
863 : Int paIndex)
864 : {
865 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
866 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
867 : CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
868 0 : Vector<Bool> deleteThem(21);
869 :
870 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
871 : Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
872 : Complex *gradVisAzData_p, *gradVisElData_p;
873 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
874 : *ConjCFMap_p, *CFMap_p;
875 : Float *l_off_p, *m_off_p;
876 : Double actualPA;
877 :
878 0 : Vector<Int> ConjCFMap, CFMap;
879 0 : actualPA = getVBPA(vb);
880 0 : ConjCFMap = polMap;
881 0 : makeCFPolMap(vb,CFMap);
882 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
883 :
884 0 : Array<Complex> rotatedConvFunc;
885 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
886 : // rotatedConvFunc,(currentCFPA-actualPA));
887 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
888 : rotatedConvFunc,0.0);
889 :
890 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
891 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
892 :
893 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
894 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
895 0 : visdata_p = visdata.getStorage(deleteThem(VISDATA));
896 0 : gradVisAzData_p = gradVisAzData.getStorage(deleteThem(GRADVISAZ));
897 0 : gradVisElData_p = gradVisElData.getStorage(deleteThem(GRADVISEL));
898 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
899 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
900 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
901 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
902 0 : dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
903 0 : vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
904 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
905 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
906 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
907 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
908 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
909 0 : vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
910 0 : vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
911 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
912 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
913 :
914 0 : Int npa=convSupport.shape()(2),actualConvSize;
915 0 : Int paIndex_Fortran = paIndex;
916 0 : actualConvSize = convFunc.shape()(0);
917 :
918 0 : IPosition shp=convSupport.shape();
919 0 : Int alwaysDoPointing=1;
920 0 : alwaysDoPointing = doPointing;
921 0 : dpbmosgrad(uvw_p,
922 : dphase_p,
923 : visdata_p,
924 0 : &s.asVector()(0),
925 0 : &s.asVector()(1),
926 : gradVisAzData_p,
927 : gradVisElData_p,
928 : &Conj,
929 : flags_p,
930 : rowFlags_p,
931 0 : &s.asVector()(2),
932 : &rownr,
933 : uvScale_p,
934 : actualOffset_p,
935 : dataPtr_p,
936 : &aNx,
937 : &aNy,
938 : &npol,
939 : &nchan,
940 : vb_freq_p,
941 : &C::c,
942 : convSupport_p,
943 : &actualConvSize,
944 : &convSampling,
945 : &wConvSize,
946 : f_convFunc_p,
947 : chanMap_p,
948 : polMap_p,
949 : &polInUse,
950 : vb_ant1_p,
951 : vb_ant2_p,
952 : &Nant_p,
953 : &scanNo,
954 : &sigma,
955 : l_off_p, m_off_p,
956 : &area,
957 : &doGrad,
958 : &alwaysDoPointing,
959 : &npa,
960 : &paIndex_Fortran,
961 : CFMap_p,
962 : ConjCFMap_p,
963 : ¤tCFPA
964 : ,&actualPA,&cfRefFreq_p
965 : );
966 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
967 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
968 :
969 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
970 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
971 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
972 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
973 0 : visdata.putStorage(visdata_p,deleteThem(VISDATA));
974 0 : gradVisAzData.putStorage(gradVisAzData_p,deleteThem(GRADVISAZ));
975 0 : gradVisElData.putStorage(gradVisElData_p,deleteThem(GRADVISEL));
976 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
977 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
978 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
979 0 : dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
980 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
981 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
982 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
983 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
984 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
985 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
986 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
987 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
988 0 : }
989 : //
990 : //----------------------------------------------------------------------
991 : //
992 0 : void PBMosaicFT::runFortranPut(Matrix<Double>& uvw,Vector<Double>& dphase,
993 : const Complex& visdata,
994 : IPosition& s,
995 : Int& /*Conj*/,
996 : Cube<Int>& flags,Vector<Int>& rowFlags,
997 : const Matrix<Float>& weight,
998 : Int& rownr,Vector<Double>& actualOffset,
999 : Array<Complex>& dataPtr,
1000 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
1001 : const VisBuffer& vb,Int& Nant_p, Int& scanNo,
1002 : Double& sigma,
1003 : Array<Float>& l_off,
1004 : Array<Float>& m_off,
1005 : Matrix<Double>& sumWeight,
1006 : Double& area,
1007 : Int& doGrad,
1008 : Int& doPSF,
1009 : Int paIndex)
1010 : {
1011 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
1012 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
1013 : CONVSUPPORT,CONVWTSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,WEIGHT,
1014 : SUMWEIGHT,CONJCFMAP,CFMAP,AVGPBPTR,CONVWTS};
1015 0 : Vector<Bool> deleteThem(25);
1016 :
1017 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
1018 : Complex *dataPtr_p, *f_convFunc_p, *f_convWts_p,*avgPBPtr;
1019 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *convWtSupport_p,
1020 : *vb_ant1_p, *vb_ant2_p,
1021 : *ConjCFMap_p, *CFMap_p;
1022 : Float *l_off_p, *m_off_p;
1023 : Float *weight_p;Double *sumwt_p,*isumwt_p;
1024 : Double actualPA;
1025 0 : const Complex *visdata_p=&visdata;
1026 :
1027 0 : Matrix<Double> iSumWt(sumWeight.shape());
1028 0 : iSumWt=0.0;
1029 0 : Vector<Int> ConjCFMap, CFMap;
1030 0 : actualPA = getVBPA(vb);
1031 :
1032 0 : ConjCFMap = polMap;
1033 0 : makeCFPolMap(vb,CFMap);
1034 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
1035 :
1036 0 : Array<Complex> rotatedConvFunc;
1037 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
1038 : // rotatedConvFunc,(currentCFPA-actualPA));
1039 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
1040 : rotatedConvFunc,0.0);
1041 :
1042 :
1043 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
1044 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
1045 :
1046 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
1047 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
1048 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
1049 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
1050 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
1051 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
1052 0 : dataPtr_p = dataPtr.getStorage(deleteThem(DATAPTR));
1053 0 : vb_freq_p = (Double *)(vb.frequency().getStorage(deleteThem(VBFREQ)));
1054 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
1055 0 : convWtSupport_p = convWtSupport.getStorage(deleteThem(CONVWTSUPPORT));
1056 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
1057 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
1058 0 : f_convWts_p = convWeights.getStorage(deleteThem(CONVWTS));
1059 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
1060 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
1061 0 : vb_ant1_p = (Int *)(vb.antenna1().getStorage(deleteThem(VBANT1)));
1062 0 : vb_ant2_p = (Int *)(vb.antenna2().getStorage(deleteThem(VBANT2)));
1063 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
1064 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
1065 0 : weight_p = (Float *)(weight.getStorage(deleteThem(WEIGHT)));
1066 0 : sumwt_p = sumWeight.getStorage(deleteThem(SUMWEIGHT));
1067 : Bool tmp;
1068 0 : isumwt_p = iSumWt.getStorage(tmp);
1069 :
1070 : // Array<Complex> avgPB_p(griddedWeights.get());
1071 0 : Array<Complex> avgPB_p;
1072 0 : if (!avgPBReady)
1073 : {
1074 0 : avgPB_p.resize(griddedWeights.shape());
1075 0 : avgPB_p=Complex(0,0);
1076 0 : avgPBPtr = avgPB_p.getStorage(deleteThem(AVGPBPTR));
1077 : }
1078 : else
1079 0 : avgPBPtr=NULL;
1080 :
1081 0 : Int npa=convSupport.shape()(2),actualConvSize, actualConvWtSize;
1082 0 : Int paIndex_Fortran = paIndex;
1083 0 : Int doAvgPB=((avgPBReady==false) &&
1084 0 : ((fabs(lastPAUsedForWtImg-actualPA)*57.2956 >= DELTAPA) ||
1085 0 : (lastPAUsedForWtImg == MAGICPAVALUE)));
1086 0 : doAvgPB=(avgPBReady==false);
1087 0 : actualConvSize = convFunc.shape()(0);
1088 0 : actualConvWtSize = convWeights.shape()(0);
1089 :
1090 0 : if (fabs(lastPAUsedForWtImg-actualPA)*57.2956 >= DELTAPA) lastPAUsedForWtImg = actualPA;
1091 :
1092 0 : IPosition shp=convSupport.shape();
1093 0 : Int alwaysDoPointing=1;
1094 0 : alwaysDoPointing = doPointing;
1095 0 : gpbmos(uvw_p,
1096 : dphase_p,
1097 : visdata_p,
1098 0 : &s.asVector()(0),
1099 0 : &s.asVector()(1),
1100 : &doPSF,
1101 : flags_p,
1102 : rowFlags_p,
1103 : weight_p,
1104 0 : &s.asVector()(2),
1105 : &rownr,
1106 : uvScale_p,
1107 : actualOffset_p,
1108 : dataPtr_p,
1109 : &aNx,
1110 : &aNy,
1111 : &npol,
1112 : &nchan,
1113 : vb_freq_p,
1114 : &C::c,
1115 : convSupport_p,
1116 : &actualConvSize,
1117 : &actualConvWtSize,
1118 : &convSampling,
1119 : &wConvSize,
1120 : f_convFunc_p,
1121 : chanMap_p,
1122 : polMap_p,
1123 : &polInUse,
1124 : // sumwt_p,
1125 : isumwt_p,
1126 : vb_ant1_p,
1127 : vb_ant2_p,
1128 : &Nant_p,
1129 : &scanNo,
1130 : &sigma,
1131 : l_off_p, m_off_p,
1132 : &area,
1133 : &doGrad,
1134 : &alwaysDoPointing,
1135 : &npa,
1136 : &paIndex_Fortran,
1137 : CFMap_p,
1138 : ConjCFMap_p,
1139 : ¤tCFPA,&actualPA,
1140 : &doAvgPB,
1141 : avgPBPtr,&cfRefFreq_p,
1142 : f_convWts_p,convWtSupport_p
1143 : );
1144 :
1145 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
1146 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
1147 :
1148 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
1149 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
1150 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
1151 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
1152 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
1153 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
1154 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
1155 0 : dataPtr.freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
1156 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
1157 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
1158 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
1159 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
1160 0 : convWeights.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVWTS));
1161 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
1162 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
1163 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
1164 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
1165 0 : weight.freeStorage((const Float*&)weight_p,deleteThem(WEIGHT));
1166 0 : sumWeight.putStorage(sumwt_p,deleteThem(SUMWEIGHT));
1167 0 : iSumWt.putStorage(isumwt_p,tmp);
1168 0 : sumWeight += iSumWt;
1169 :
1170 0 : if (!avgPBReady)
1171 : {
1172 0 : nApertures+=Complex(1.0,0.0);
1173 : // Get the griddedWeigths as a referenced array
1174 0 : Array<Complex> gwts; Bool removeDegenerateAxis=false;
1175 0 : griddedWeights.get(gwts, removeDegenerateAxis);
1176 : // griddedWeights.put(griddedWeights.get()+avgPB_p);
1177 0 : gwts = gwts + avgPB_p;
1178 : // if (!reference)
1179 : // griddedWeights.put(gwts);
1180 : }
1181 0 : }
1182 : //
1183 : //---------------------------------------------------------------
1184 : //
1185 : // Initialize the FFT to the Sky. Here we have to setup and
1186 : // initialize the grid.
1187 : //
1188 0 : void PBMosaicFT::initializeToSky(ImageInterface<Complex>& iimage,
1189 : Matrix<Float>& weight,
1190 : const VisBuffer& vb)
1191 : {
1192 0 : image=&iimage;
1193 :
1194 0 : init();
1195 0 : initMaps(vb);
1196 0 : nx = image->shape()(0);
1197 0 : ny = image->shape()(1);
1198 0 : npol = image->shape()(2);
1199 0 : nchan = image->shape()(3);
1200 :
1201 0 : if(image->shape().product()>cachesize) isTiled=true;
1202 0 : else isTiled=false;
1203 :
1204 :
1205 0 : sumWeight=0.0;
1206 0 : weight.resize(sumWeight.shape());
1207 0 : weight=0.0;
1208 :
1209 0 : if(isTiled)
1210 : {
1211 0 : imageCache->flush();
1212 0 : image->set(Complex(0.0));
1213 : //lattice=image;
1214 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
1215 : }
1216 : else
1217 : {
1218 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1219 0 : griddedData.resize(gridShape);
1220 0 : griddedData=Complex(0.0);
1221 : // if(arrayLattice) delete arrayLattice; arrayLattice=0;
1222 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1223 0 : lattice=arrayLattice;
1224 : }
1225 : //AlwaysAssert(lattice, AipsError);
1226 0 : PAIndex = -1;
1227 0 : if (resetPBs)
1228 : {
1229 0 : griddedWeights.resize(iimage.shape());
1230 0 : griddedWeights.setCoordinateInfo(iimage.coordinates());
1231 0 : griddedWeights.set(0.0);
1232 0 : noOfPASteps = 1;
1233 0 : pbPeaks.resize(griddedWeights.shape()(2));
1234 0 : pbPeaks.set(0.0);
1235 0 : resetPBs=false;
1236 : }
1237 0 : }
1238 : //
1239 : //---------------------------------------------------------------
1240 : //
1241 0 : void PBMosaicFT::finalizeToSky()
1242 : {
1243 0 : if(isTiled)
1244 : {
1245 0 : logIO() << LogOrigin("PBMosaicFT", "finalizeToSky") << LogIO::NORMAL;
1246 :
1247 0 : AlwaysAssert(image, AipsError);
1248 0 : AlwaysAssert(imageCache, AipsError);
1249 0 : imageCache->flush();
1250 0 : ostringstream o;
1251 0 : imageCache->showCacheStatistics(o);
1252 0 : logIO() << o.str() << LogIO::POST;
1253 : }
1254 :
1255 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
1256 0 : PAIndex = -1;
1257 :
1258 0 : paChangeDetector.reset();
1259 0 : cfCache.finalize();
1260 0 : convFuncCacheReady=true;
1261 0 : }
1262 :
1263 : } //# NAMESPACE CASA - END
|
