Line data Source code
1 : // -*- C++ -*-
2 : //# CFBuffer.cc: Implementation of the CFBuffer 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 : #include <synthesis/TransformMachines/CFBuffer.h>
29 : #include <synthesis/TransformMachines/Utils.h>
30 : #include <casacore/casa/Utilities/BinarySearch.h>
31 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
32 : #include <casacore/casa/Utilities/Assert.h>
33 : using namespace casacore;
34 :
35 : namespace casa{
36 :
37 : // Instantiate a commonly used extern template
38 :
39 : //
40 : //---------------------------------------------------------------
41 : //
42 : //
43 0 : void CFBuffer::setParams(const CFBuffer& other)
44 : {
45 0 : wValues_p.assign(other.wValues_p);
46 0 : freqValues_p.assign(other.freqValues_p);
47 0 : muellerElements_p.assign(other.muellerElements_p);
48 0 : muellerElementsIndex_p.assign(other.muellerElementsIndex_p);
49 0 : conjMuellerElements_p.assign(other.conjMuellerElements_p);
50 0 : conjMuellerElementsIndex_p.assign(other.conjMuellerElementsIndex_p);
51 0 : wValIncr_p = other.wValIncr_p;
52 0 : freqValIncr_p = other.freqValIncr_p;
53 0 : muellerMask_p.assign(other.muellerMask_p);
54 0 : pointingOffset_p.assign(other.pointingOffset_p);
55 0 : freqNdxMapsReady_p = other.freqNdxMapsReady_p;
56 :
57 0 : freqNdxMap_p.assign(other.freqNdxMap_p);
58 0 : for(uInt i=0;i<freqNdxMap_p.nelements();i++) freqNdxMap_p[i].assign(other.freqNdxMap_p[i]);
59 0 : conjFreqNdxMap_p.assign(other.conjFreqNdxMap_p);
60 0 : for(uInt i=0;i<conjFreqNdxMap_p.nelements();i++) conjFreqNdxMap_p[i].assign(other.conjFreqNdxMap_p[i]);
61 0 : cfCacheDirName_p = other.cfCacheDirName_p;
62 0 : }
63 : //---------------------------------------------------------------
64 : //
65 : //
66 0 : CountedPtr<CFBuffer> CFBuffer::clone()
67 : {
68 0 : CountedPtr<CFBuffer> clone=new CFBuffer();
69 0 : clone->setParams(*this);
70 :
71 0 : IPosition shp(cfCells_p.shape());
72 0 : clone->resize(shp);
73 0 : clone->allocCells(cfCells_p);
74 : // clone->show("####CLONE: ",cerr);
75 0 : return clone;
76 : }
77 0 : void CFBuffer::allocCells(const Cube<CountedPtr<CFCell> >& cells)
78 : {
79 0 : IPosition shp(cells.shape());
80 0 : for (Int i=0;i < shp[0]; i++)
81 0 : for (Int j=0; j < shp[1]; j++)
82 0 : for (Int k=0; k < shp[2]; k++)
83 : {
84 0 : cfCells_p(i,j,k) = cells(i,j,k)->clone();
85 : }
86 0 : }
87 : //---------------------------------------------------------------
88 : //
89 : // template<class T> Int CFBuffer<T>
90 0 : Int CFBuffer::noOfMuellerElements(const PolMapType& muellerElements)
91 : {
92 0 : Int n=0,nrows = muellerElements.nelements();
93 0 : for (Int i=0;i<nrows;i++)
94 0 : n+=muellerElements(i).nelements();
95 0 : return n;
96 : // Int n=0;
97 : // for(Int i=0;i<muellerElements.nelements();i++)
98 : // for(Int j=0;j<muellerElements(i).nelements();j++)
99 : // if (muellerElements(i)(j) > n) n=muellerElements(i)(j);
100 : // return n;
101 : }
102 : //
103 : //---------------------------------------------------------------
104 : //
105 0 : void CFBuffer::clear()
106 : {
107 0 : IPosition shp(cfCells_p.shape());
108 0 : for (Int i=0;i < shp[0]; i++)
109 0 : for (Int j=0; j < shp[1]; j++)
110 0 : for (Int k=0; k < shp[2]; k++)
111 : {
112 0 : cfCells_p(i,j,k)->clear();
113 : }
114 0 : }
115 : //
116 : //---------------------------------------------------------------
117 : //
118 : // template<class T> void CFBuffer<T>
119 0 : void CFBuffer::resize(const Double& wScale, const Double& freqIncr,
120 : const Vector<Double>& wValues,
121 : const Vector<Double>& freqValues,
122 : const PolMapType& muellerElements,
123 : const PolMapType& muellerElementsIndex,
124 : const PolMapType& conjMuellerElements,
125 : const PolMapType& conjMuellerElementsIndex
126 : )
127 : {
128 0 : wValues_p.assign(wValues);
129 0 : freqValues_p.assign(freqValues);
130 0 : wValIncr_p = wScale;
131 0 : freqValIncr_p = freqIncr;
132 : // muellerMask_p.assign(muellerElements);
133 :
134 : // nPol_p=noOfMuellerElements(muellerMask_p);
135 0 : nPol_p=noOfMuellerElements(muellerElementsIndex);
136 0 : nChan_p = freqValues_p.nelements();
137 0 : nW_p = wValues_p.nelements();
138 :
139 : // muellerElements_p.resize(nPol_p);
140 0 : muellerElements_p.assign(muellerElements);
141 0 : muellerElementsIndex_p.assign(muellerElementsIndex);
142 0 : conjMuellerElements_p.assign(conjMuellerElements);
143 0 : conjMuellerElementsIndex_p.assign(conjMuellerElementsIndex);
144 : // Resize the various aux. information storage buffers.
145 : //
146 : // Resize the storage. Retain the value of the existing pixels.
147 : // New pixels due to resize, if any, are assigned a new Array<T>
148 : // pointer.
149 0 : cfCells_p.resize(nChan_p, nW_p, nPol_p, true);
150 :
151 0 : for (uInt i=0;i<cfCells_p.shape()(0);i++) // nChan_p
152 0 : for (uInt j=0;j<cfCells_p.shape()(1);j++) // nW_p
153 : {
154 0 : Int k=0; // nPol_p
155 0 : for(uInt prow=0;prow<muellerElements_p.nelements();prow++)
156 0 : for(uInt pcol=0;pcol<muellerElements_p(prow).nelements();pcol++)
157 : {
158 0 : if (cfCells_p(i,j,k).null()) cfCells_p(i,j,k) = new CFCell;
159 0 : if (cfCells_p(i,j,k)->storage_p.null()) cfCells_p(i,j,k)->storage_p=new Array<TT>;
160 0 : cfCells_p(i,j,k)->freqValue_p = freqValues(i);
161 0 : cfCells_p(i,j,k)->freqIncr_p = freqIncr;
162 0 : cfCells_p(i,j,k)->wValue_p = wValues(j);
163 0 : cfCells_p(i,j,k)->wIncr_p = wValIncr_p;
164 0 : cfCells_p(i,j,k)->muellerElement_p = muellerElements_p(prow)(pcol);
165 0 : cfCells_p(i,j,k)->xSupport_p = 0;
166 0 : cfCells_p(i,j,k)->ySupport_p = 0;
167 0 : k++;
168 : }
169 : }
170 0 : }
171 : //
172 : //---------------------------------------------------------------
173 : //
174 : // template <class T> void CFBuffer<T>
175 : // RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& muellerElement,
176 : // const TableRecord& miscInfo)
177 : // {
178 : // RigidVector<Int,3> ndx = setParams(inu, iw,
179 : // miscInfo.coordSys, miscInfo.sampling,
180 : // miscInfo.xSupport,miscInfo.ySupport,
181 : // miscInfo.freqValue, miscInfo.wValue,
182 : // muellerElements,
183 : // miscInfo.fileName, miscInfo.conjfFreq,
184 : // miscinfo.conjPoln);
185 : // cfCells_p(ndx(0), ndx(1), ndx(2))->telescopeName = micsInfo.telescopeName;
186 : // }
187 : //
188 : //---------------------------------------------------------------
189 : //
190 : // template <class T> void CFBuffer<T>
191 0 : RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& ipx, const Int& ipy,
192 : const Double& freqValue,
193 : const Double& wValue,
194 : const Int& muellerElement,
195 : CoordinateSystem& cs,
196 : const TableRecord& miscInfo)
197 : {
198 0 : float sampling; miscInfo.get("Sampling",sampling);
199 0 : int xSupport, ySupport; miscInfo.get("Xsupport",xSupport);miscInfo.get("Ysupport",ySupport);
200 0 : String fileName; miscInfo.get("Name",fileName);
201 0 : double conjFreq; miscInfo.get("ConjFreq", conjFreq);
202 0 : int conjPoln; miscInfo.get("ConjPoln", conjPoln);
203 0 : String telescopeName; miscInfo.get("TelescopeName", telescopeName);
204 0 : String bandName; miscInfo.get("BandName", bandName);
205 0 : float diameter; miscInfo.get("Diameter", diameter);
206 : // In the absense of evidence, assume that users are sensible and
207 : // are using AWProjection where it is really need it and not for
208 : // using it as a replacement for rotatially symmetric stuff. So
209 : // by default, the CFs are assumed to be rotationally asymmetric.
210 0 : bool isRotationallySymmetric=True;
211 :
212 0 : if (miscInfo.isDefined("OpCode"))
213 0 : miscInfo.get("OpCode",isRotationallySymmetric);
214 :
215 : RigidVector<Int,3> ndx=setParams(inu, iw, ipx, ipy, freqValue, bandName, wValue, muellerElement, cs,
216 : sampling, xSupport, ySupport, fileName, conjFreq, conjPoln, telescopeName,
217 0 : diameter);
218 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->isRotationallySymmetric_p = isRotationallySymmetric;
219 0 : return ndx;
220 : }
221 0 : RigidVector<Int, 3> CFBuffer::setParams(const Int& inu, const Int& iw, const Int& /*ipx*/, const Int& /*ipy*/,
222 : const Double& freqValue,
223 : const String& bandName,
224 : const Double& wValue,
225 : const Int& muellerElement,
226 : CoordinateSystem& cs,
227 : Float& sampling,
228 : Int& xSupport, Int& ySupport,
229 : const String& fileName,
230 : const Double& conjFreq,
231 : const Int& conjPoln,
232 : const String& telescopeName,
233 : const Float& diameter)
234 : {
235 0 : RigidVector<Int,3> ndx=getIndex(freqValue, wValue, muellerElement);
236 0 : ndx(0)=inu; ndx(1)=iw;//ndx(2) = muellerElements_p(ipx)(ipy);
237 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->sampling_p = sampling;
238 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->xSupport_p = xSupport;
239 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->ySupport_p = ySupport;
240 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->muellerElement_p = muellerElement;
241 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->conjFreq_p = conjFreq;
242 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->conjPoln_p = conjPoln;
243 0 : if (fileName != "")
244 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->fileName_p = fileName;
245 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->telescopeName_p = telescopeName;//String("EVLA");
246 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->diameter_p = diameter;//25.0;
247 0 : if (bandName != "") cfCells_p(ndx(0),ndx(1),ndx(2))->bandName_p = bandName;
248 :
249 0 : Int index=cs.findCoordinate(Coordinate::SPECTRAL);
250 0 : SpectralCoordinate spCS = cs.spectralCoordinate(index);
251 0 : Vector<Double> val; val.resize(1);val(0)=freqValues_p(ndx(0));
252 0 : spCS.setReferenceValue(val);
253 : //val.resize();
254 : //val = spCS.increment();
255 : //val = cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p;
256 : //spCS.setIncrement(val);
257 0 : cs.replaceCoordinate(spCS,index);
258 :
259 : // cfCells_p(ndx(0),ndx(1),ndx(2))->freqValue_p =
260 : // cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)).referenceValue()(0);
261 : // cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p =
262 : // cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)).increment()(0);
263 :
264 0 : cfCells_p(ndx(0),ndx(1),ndx(2))->coordSys_p = cs;
265 0 : return ndx;
266 : }
267 : //
268 : //---------------------------------------------------------------
269 : //
270 0 : void CFBuffer::setPA(Float& pa)
271 : {
272 0 : RigidVector<Int,3> ndx(-1);
273 0 : Int nF=cfCells_p.shape()(0), nW=cfCells_p.shape()(1), nM=cfCells_p.shape()(2);
274 0 : for (Int i=0; i<nF; i++)
275 0 : for(Int j=0; j<nW; j++)
276 0 : for(Int k=0; k<nM; k++)
277 0 : cfCells_p(i,j,k)->pa_p=Quantity(pa,"deg");
278 0 : }
279 :
280 0 : void CFBuffer::setParams(Int& /*nx*/, Int& /*ny*/, CoordinateSystem& /*cs*/, Float& /*sampling*/,
281 : Int& /*xSupport*/, Int& /*ySupport*/, const Double& /*freqValue*/,
282 : const Double& /*wValue*/, const Int& /*muellerElement*/,
283 : const String& /*fileName*/)
284 : {
285 : /*
286 : RigidVector<Int,3> ndx=setParams(cs, sampling, xSupport, ySupport,
287 : freqValue, wValue, muellerElement);
288 : // Adding degenerate axis here to make life easier when using these
289 : // buffers to write them as Images later (with a 4D co-ordinate
290 : // system - the degenerate axis become the Freq. and Polarization
291 : // axis).
292 :
293 : // cfCells_p(ndx(0),ndx(1),ndx(2))->storage_p->resize(IPosition(4,nx,ny,1,1));
294 : cfCells_p(ndx(0),ndx(1),ndx(2))->shape_p=IPosition(4,nx,ny,1,1);
295 : */
296 0 : }
297 : //
298 : //---------------------------------------------------------------
299 : //
300 : // template<class T> Array<T>& CFBuffer<T>::
301 0 : CFCell& CFBuffer::getCFCell(const Int& i, const Int& j, const Int& k)
302 : {
303 0 : return *cfCells_p(i,j,k); // Nfreq x Nw x Nmueller
304 : }
305 : //
306 : //---------------------------------------------------------------
307 : //
308 : // template <class T> Array<T>& CFBuffer<T>
309 0 : CFCell& CFBuffer::getCFCell(const Double& freqVal,
310 : const Double& wValue,
311 : const Int& muellerElement) // muellerElement: (i,j) of the Mueller Matrix
312 : {
313 0 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
314 0 : return getCFCell(ndx(0), ndx(1),ndx(2));
315 : }
316 : //
317 : //---------------------------------------------------------------
318 : //
319 : // template<class T> Array<T>& CFBuffer<T>::
320 0 : CountedPtr<CFCell>& CFBuffer::getCFCellPtr(const Int& i/*FreqNdx*/, const Int& j/*WNdx*/, const Int& k/*MuellerNdx*/)
321 : {
322 : // IPosition shp=cfCells_p.shape();
323 : // if (cfHitsStats.shape().product()==0)
324 : // {
325 : // cfHitsStats.resize(shp);
326 : // cfHitsStats = 0;
327 : // }
328 : // try
329 : // {
330 : // AlwaysAssert(((i<shp[0]) && (j<shp[1]) && (k<shp[2])) , AipsError);
331 : // }
332 : // catch (AipsError x)
333 : // {
334 : // cerr << "#### " << i << " " << j << " " << k << " " << shp << endl;
335 : // throw(x);
336 : // }
337 :
338 : // cfHitsStats(i,j,k)++;
339 : // // cerr << "CFBuffer: Cell: " << i << " " << j << " " << k << " " << cfCells_p(i,j,k)->xSupport_p << endl;
340 0 : return cfCells_p(i,j,k); // Nfreq x Nw x Nmueller
341 : }
342 : //
343 : //---------------------------------------------------------------
344 : //
345 : // template <class T> Array<T>& CFBuffer<T>
346 0 : CountedPtr<CFCell>& CFBuffer::getCFCellPtr(const Double& freqVal,
347 : const Double& wValue,
348 : const Int & muellerElement)
349 : {
350 0 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
351 0 : return cfCells_p(ndx(0), ndx(1),ndx(2));
352 : }
353 : //
354 : //---------------------------------------------------------------
355 : //
356 : // template<class T> Vector<Int>& CFBuffer<T>
357 : // This should be made more efficient with binary search.
358 0 : RigidVector<Int,3> CFBuffer::getIndex(const Double& freqValue,
359 : const Double& wValue,
360 : const Int& muellerElement)
361 : {
362 0 : RigidVector<Int,3> ndx(-1);
363 0 : Int nF=cfCells_p.shape()(0), nW=cfCells_p.shape()(1), nM=cfCells_p.shape()(2);
364 : Int i,j,k;
365 : //UNUSED: Int di;
366 :
367 : // Double dfMin=0;di=0;
368 : // for (i=0;i<nF;i++)
369 : // {
370 : // Double df=fabs(cfCells_p(i,0,0)->freqValue_p - freqValue);
371 : // if (df < dfMin) {dfMin=df; di=i;}
372 : // }
373 : // i=di;
374 0 : for (i=0; i<nF; i++)
375 : // if ((fabs(cfCells_p(i,0,0)->freqValue_p - freqValue) < cfCells_p(i,0,0)->freqIncr_p))
376 0 : if (cfCells_p(i,0,0)->freqValue_p == freqValue) break;
377 0 : for(j=0; j<nW; j++) if (cfCells_p(i,j,0)->wValue_p == wValue) break;
378 0 : for (k=0; k<nM; k++) if (cfCells_p(i,j,k)->muellerElement_p == muellerElement) break;
379 :
380 0 : ndx(0)=i;ndx(1)=j;ndx(2)=k;
381 : // for (Int i=0;i<nF;i++)
382 : // for(Int j=0;j<nW;j++)
383 : // for(Int k=0;k<nM;k++)
384 : // {
385 : // if (
386 : // (fabs(cfCells_p(i,j,k)->freqValue_p - freqValue) < cfCells_p(i,j,k)->freqIncr_p) &&
387 : // (cfCells_p(i,j,k)->wValue_p == wValue) &&
388 : // (cfCells_p(i,j,k)->muellerElement_p == muellerElement)
389 : // )
390 : // {ndx(0)=i;ndx(1)=j;ndx(2)=k;break;}
391 : // }
392 :
393 : // cerr << "@#$%#@ " << freqValue
394 : // << " " << cfCells_p(ndx(0),ndx(1),ndx(2))->freqValue_p
395 : // << " " << cfCells_p(ndx(0),ndx(1),ndx(2))->freqIncr_p
396 : // << " " << ndx
397 : // << endl;
398 : // for(uInt i=0;i<nF;i++)
399 : // if (freqValue==cfCells_pfreqValues_p(i)) {ndx(0)=i;break;}
400 :
401 : // for(uInt i=0;i<wValues_p.nelements();i++)
402 : // if (wValue==wValues_p(i)) {ndx(1)=i;break;}
403 :
404 : // ndx(2)=muellerElements_p(muellerElement(0),muellerElement(1));
405 :
406 0 : return ndx;
407 : }
408 : //
409 : //---------------------------------------------------------------
410 : //
411 0 : void CFBuffer::getParams(CoordinateSystem& cs, Float& sampling,
412 : Int& xSupport, Int& ySupport,
413 : const Double& freqVal, const Double& wValue,
414 : const Int& muellerElement)
415 : {
416 0 : RigidVector<Int,3> ndx=getIndex(freqVal, wValue, muellerElement);
417 0 : getParams(cs, sampling, xSupport, ySupport, ndx(0), ndx(1), ndx(2));
418 0 : }
419 : //
420 : //---------------------------------------------------------------
421 : //
422 0 : Double CFBuffer::nearest(Bool& found, const Double& val,
423 : const Vector<Double>& valList,
424 : const Double& /*incr*/)
425 : {
426 0 : Int n=valList.nelements();
427 0 : if (n==1) {found=true;return valList[0];}
428 0 : Int where = binarySearch(found, valList, val, n);
429 0 : if (found) return valList(where);
430 0 : else return -1.0;
431 : }
432 : //
433 : //---------------------------------------------------------------
434 : //
435 0 : Int CFBuffer::nearestNdx(const Double& val,
436 : const Vector<Double>& /*valList*/,
437 : const Double& incr)
438 : {
439 0 : return SynthesisUtils::nint(incr*val);
440 :
441 : // Int n=valList.nelements();
442 : // if (n==1) return 0;
443 : // else return -1;
444 :
445 : // Int where = binarySearch(found, valList, val, n);
446 : // if (found) return valList(where);
447 : // else return -1.0;
448 : }
449 : //
450 : //---------------------------------------------------------------
451 : //
452 : // template <class T> void CFBuffer<T>
453 0 : void CFBuffer::show(const char *Mesg,ostream &os)
454 : {
455 0 : LogIO log_l(LogOrigin("CFBuffer","show[R&D]"));
456 :
457 0 : if (Mesg != NULL) os << Mesg << endl;
458 0 : os<< "---------------------------------------------------------" << endl
459 0 : << "Shapes: " << cfCells_p.shape() << endl;
460 0 : for (Int i=0;i<cfCells_p.shape()(0);i++)
461 0 : for (Int j=0;j<cfCells_p.shape()(1);j++)
462 0 : for (Int k=0;k<cfCells_p.shape()(2);k++)
463 : {
464 0 : os << "CFCell["<<i<<","<<j<<","<<k<<"]:" << endl;
465 0 : cfCells_p(i,j,k)->show(Mesg, os);
466 0 : os << "Pointing offset: " << pointingOffset_p << endl;
467 : }
468 0 : os << "---------------------------------------------------------" << endl;
469 0 : }
470 :
471 0 : void CFBuffer::makePersistent(const char *dir, const char *cfName)
472 : {
473 0 : for (Int i=0;i<cfCells_p.shape()(0);i++)
474 0 : for (Int j=0;j<cfCells_p.shape()(1);j++)
475 0 : for (Int k=0;k<cfCells_p.shape()(2);k++)
476 : {
477 0 : ostringstream name;
478 0 : name << String(cfName) << "_CF_" << i << "_" << j << "_" << k << ".im";
479 : //cerr << "CFB Name : " << name.str() << endl;
480 : // const char *formedName;
481 : // if (cfName != "" ) formedName = name.str().c_str();
482 : // else formedName = cfName;
483 : // cerr << "Formed name = " << formedName << endl;
484 0 : cfCells_p(i,j,k)->makePersistent(dir, name.str().c_str());
485 : }
486 0 : }
487 :
488 0 : Int CFBuffer::nearestFreqNdx(const Double& freqVal)
489 : {
490 : Int index;
491 0 : SynthesisUtils::nearestValue(freqValues_p, freqVal, index);
492 0 : return index;
493 : // // The algorithm below has a N*log(N) cost.
494 : // Vector<Double> diff = fabs(freqValues_p - freqVal);
495 : // Bool dummy;
496 : // Sort sorter(diff.getStorage(dummy), sizeof(Double));
497 : // sorter.sortKey((uInt)0,TpDouble);
498 : // Int nch=freqValues_p.nelements();
499 : // Vector<uInt> sortIndx;
500 : // sorter.sort(sortIndx, nch);
501 :
502 : // return sortIndx(0);
503 :
504 : // Int ndx=min(freqValues_p.nelements()-1,max(0,SynthesisUtils::nint((freqVal-freqValues_p[0])/freqValIncr_p)));
505 : // return ndx;
506 : }
507 :
508 :
509 0 : void CFBuffer::primeTheCache()
510 : {
511 : //
512 : // In each CFCell of this CFBuffer, cache things that might be
513 : // required in tight loops and can be expensive to extract
514 : // otherwise.
515 : //
516 0 : IPosition cfCShape=cfCells_p.shape();
517 0 : for (Int i=0; i < cfCShape(0); i++)
518 0 : for (Int j=0; j < cfCShape(1); j++)
519 0 : for (Int k=0; k < cfCShape(2); k++)
520 0 : getCFCellPtr(i,j,k)->initCache();
521 0 : }
522 :
523 0 : void CFBuffer::initMaps(const VisBuffer&, // vb,
524 : const Matrix<Double>& freqSelection, const Double& imRefFreq)
525 : {
526 0 : Vector<Double> spwList=freqSelection.column(0);
527 0 : Int maxSpw=(Int)(max(spwList));
528 0 : freqNdxMap_p.resize(maxSpw+1);
529 0 : conjFreqNdxMap_p.resize(maxSpw+1);
530 :
531 0 : for (Int i=0;i<(Int)spwList.nelements(); i++)
532 : {
533 0 : Int spw=(Int)freqSelection(i,0);
534 0 : Double fmin=freqSelection(i,1), fmax=freqSelection(i,2), finc=freqSelection(i,3);
535 0 : Int nchan = (Int)((fmax-fmin)/finc + 1);
536 0 : freqNdxMap_p[spw].resize(nchan);
537 0 : conjFreqNdxMap_p[spw].resize(nchan);
538 0 : for (Int c=0;c<nchan;c++)
539 : {
540 0 : Double freq=fmin+c*finc;
541 0 : Double conjFreq=sqrt(2*imRefFreq*imRefFreq - freq*freq);
542 0 : freqNdxMap_p[spw][c]=nearestFreqNdx(freq);
543 0 : conjFreqNdxMap_p[spw][c]=nearestFreqNdx(conjFreq);
544 : }
545 : }
546 :
547 :
548 : // cerr << "CFBuffer::initMaps: "
549 : // << freqSelection << endl
550 : // << freqValues_p << endl
551 : // << freqNdxMap_p << endl
552 : // << conjFreqNdxMap_p << endl;
553 :
554 0 : }
555 :
556 0 : void CFBuffer::initPolMaps(PolMapType& polMap, PolMapType& conjPolMap)
557 : {
558 0 : muellerElementsIndex_p.assign(polMap);
559 0 : conjMuellerElementsIndex_p.assign(conjPolMap);
560 0 : }
561 :
562 0 : void CFBuffer::getFreqNdxMaps(Vector<Vector<Int> >& freqNdx, Vector<Vector<Int> >& conjFreqNdx)
563 : {
564 : Int nspw;
565 0 : nspw=freqNdxMap_p.nelements();
566 0 : freqNdx.resize(nspw);
567 0 : for (Int s=0;s<nspw;s++)
568 0 : freqNdx[s].assign(freqNdxMap_p[s]);
569 :
570 0 : nspw=conjFreqNdxMap_p.nelements();
571 0 : conjFreqNdx.resize(nspw);
572 0 : for (Int s=0;s<nspw;s++)
573 0 : conjFreqNdx[s].assign(conjFreqNdxMap_p[s]);
574 0 : }
575 :
576 0 : void CFBuffer::ASSIGNVVofI(Int** &target,Vector<Vector<Int> >& source, Bool& doAlloc)
577 : {
578 : // cerr << "Source: " << target << " " << source << endl;
579 :
580 : Int nx,ny;
581 0 : Bool b=(doAlloc||(target==NULL));
582 0 : nx=source.nelements();
583 0 : if (b) target=(int **)malloc(nx*sizeof(int*));
584 0 : for (int i=0;i<nx;i++)
585 : {
586 0 : ny = source[i].nelements();
587 0 : if (b) (target[i]) = (int *)malloc(ny*sizeof(int));
588 0 : for (int j=0;j<ny;j++)
589 0 : (target)[i][j]=source[i][j];
590 : }
591 :
592 : // cerr << "Target: ";
593 : // for (int ii=0;ii<source.nelements();ii++)
594 : // {
595 : // Int ny=source[ii].nelements();
596 : // for(Int jj=0;jj<ny;jj++)
597 : // cerr << target[ii][jj] << " ";
598 : // cerr << endl;
599 : // }
600 0 : }
601 :
602 0 : void CFBuffer::getAsStruct(CFBStruct& st)
603 : {
604 0 : Vector<Int> shp=cfCells_p.shape().asVector();
605 :
606 0 : Bool doAlloc=(st.CFBStorage == NULL);
607 :
608 0 : st.shape[0]=shp[0];
609 0 : st.shape[1]=shp[1];
610 0 : st.shape[2]=shp[2];
611 0 : st.fIncr = freqValIncr_p; st.wIncr = wValIncr_p;
612 :
613 :
614 : Bool dummy;
615 0 : Int n=cfCells_p.nelements();
616 0 : if (doAlloc) st.CFBStorage = (CFCStruct *)malloc(n*sizeof(CFCStruct));
617 0 : CountedPtr<CFCell> *cfstore=cfCells_p.getStorage(dummy);
618 0 : for (int i=0;i<n;i++)
619 : {
620 : // st.CFBStorage[i]=(CFCStruct *)malloc(sizeof(CFCStruct));
621 0 : (cfstore[i]).operator->()->getAsStruct(st.CFBStorage[i]);
622 : }
623 : // st.CFBStorage[i] = (cfstore[i]).operator->()->getStorage()->getStorage(dummy);
624 :
625 0 : if (doAlloc) st.pointingOffset=(Double *)malloc(pointingOffset_p.nelements()*sizeof(Double));
626 0 : for (uInt i=0;i<pointingOffset_p.nelements();i++) st.pointingOffset[i]=pointingOffset_p[i];
627 :
628 0 : if (doAlloc) st.freqValues=(Double *)malloc(freqValues_p.nelements()*sizeof(Double));
629 0 : for (uInt i=0;i<freqValues_p.nelements();i++) st.freqValues[i]=freqValues_p[i];
630 :
631 0 : if (doAlloc) st.wValues=(Double *)malloc(wValues_p.nelements()*sizeof(Double));
632 0 : for (uInt i=0;i<wValues_p.nelements();i++) st.wValues[i]=wValues_p[i];
633 :
634 0 : ASSIGNVVofI(st.muellerElements, muellerElements_p, doAlloc);
635 0 : ASSIGNVVofI(st.muellerElementsIndex, muellerElementsIndex_p,doAlloc);
636 0 : ASSIGNVVofI(st.conjMuellerElements, conjMuellerElements_p,doAlloc);
637 0 : ASSIGNVVofI(st.conjMuellerElementsIndex, conjMuellerElementsIndex_p,doAlloc);
638 :
639 0 : st.nMueller = muellerElements_p.nelements();
640 :
641 0 : Vector<Vector<Int> > freqNdx, conjFreqNdx;
642 0 : getFreqNdxMaps(freqNdx, conjFreqNdx);
643 0 : ASSIGNVVofI(st.freqNdxMap, freqNdx, doAlloc);
644 0 : ASSIGNVVofI(st.conjFreqNdxMap, conjFreqNdx, doAlloc);
645 0 : }
646 :
647 : //
648 : //----------------------------------------------------------------------
649 : //
650 0 : void CFBuffer::fill(const Int& /*nx*/, const Int& /*ny*/,
651 : const Vector<Double>& freqValues,
652 : const Vector<Double>& wValues,
653 : const PolMapType& muellerElements)
654 : {
655 :
656 0 : LogIO log_l(LogOrigin("CFBuffer", "fillBuffer[R&D]"));
657 0 : for (uInt imx=0;imx<muellerElements.nelements();imx++) // Loop over all MuellerElements
658 0 : for (uInt imy=0;imy<muellerElements(imx).nelements();imy++)
659 : {
660 0 : for (uInt inu=0;inu<freqValues.nelements();inu++) // All freq. channels
661 : {
662 0 : for (uInt iw=0;iw<wValues.nelements();iw++) // All w-planes
663 : {
664 : log_l << " CF("
665 0 : << "M:"<<muellerElements(imx)(imy)
666 : << ",C:" << inu
667 0 : << ",W:" << iw << "): ";
668 0 : Vector<Double> ftRef(2);
669 :
670 : // ftRef(0)=cfWtBuf.shape()(0)/2.0;
671 : // ftRef(1)=cfWtBuf.shape()(1)/2.0;
672 0 : CoordinateSystem ftCoords;//=cs_l;
673 : // SynthesisUtils::makeFTCoordSys(cs_l, cfWtBuf.shape()(0), ftRef, ftCoords);
674 :
675 : // cfb.setParams(inu,iw,imx,imy,//muellerElements(imx)(imy),
676 : // ftCoords, sampling, xSupport, ySupport,
677 : // freqValues(inu), wValues(iw), muellerElements(imx)(imy));
678 : // cfb.getCFCellPtr(freqValues(inu), wValues(iw),
679 : // muellerElements(imx)(imy))->pa_p=Quantity(vbPA,"rad");
680 : //
681 : // Now tha the CFs have been computed, cache its
682 : // paramters in CFCell for quick access in tight
683 : // loops (in the re-sampler, e.g.).
684 : //
685 0 : (getCFCellPtr(freqValues(inu), wValues(iw),
686 0 : muellerElements(imx)(imy)))->initCache();
687 : }
688 : }
689 : }
690 0 : }
691 :
692 : } // end casa namespace
693 :
694 :
695 :
|
