Line data Source code
1 : // -*- C++ -*-
2 : //# Utils.cc: Implementation of global functions from Utils.h
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 <casacore/ms/MeasurementSets/MSRange.h>
29 : #include <msvis/MSVis/VisBuffer2.h>
30 : #include <casacore/casa/Logging/LogIO.h>
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <casacore/measures/Measures/MEpoch.h>
33 : #include <casacore/measures/Measures/MeasTable.h>
34 : #include <synthesis/TransformMachines2/Utils.h>
35 : #include <synthesis/TransformMachines/StokesImageUtil.h>
36 : #include <casacore/casa/Utilities/Assert.h>
37 : #include <casacore/casa/Arrays/Vector.h>
38 : #include <casacore/casa/Arrays/ArrayMath.h>
39 : #include <casacore/lattices/LEL/LatticeExpr.h>
40 : #include <casacore/images/Images/PagedImage.h>
41 : #include <casacore/images/Images/ImageRegrid.h>
42 : #include <casacore/casa/Containers/Record.h>
43 : #include <casacore/lattices/Lattices/LatticeIterator.h>
44 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
45 : #include <casacore/lattices/Lattices/LatticeStepper.h>
46 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
47 : #include <casacore/casa/System/Aipsrc.h>
48 : #include <msvis/MSVis/VisibilityIterator2.h>
49 :
50 : using namespace casacore;
51 : namespace casa{
52 : using namespace vi;
53 : namespace refim{
54 : //
55 : //--------------------------------------------------------------------------------------------
56 : //
57 0 : void storeImg(String fileName,ImageInterface<Complex>& theImg, Bool writeReIm)
58 : {
59 0 : PagedImage<Complex> ctmp(theImg.shape(), theImg.coordinates(), fileName);
60 0 : LatticeExpr<Complex> le(theImg);
61 0 : ctmp.copyData(le);
62 0 : if (writeReIm)
63 : {
64 0 : ostringstream reName,imName;
65 0 : reName << "re" << fileName;
66 0 : imName << "im" << fileName;
67 : {
68 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), reName);
69 : //LatticeExpr<Float> le(abs(theImg));
70 0 : LatticeExpr<Float> le(real(theImg));
71 0 : tmp.copyData(le);
72 : }
73 : {
74 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), imName);
75 0 : LatticeExpr<Float> le(arg(theImg));
76 0 : tmp.copyData(le);
77 : }
78 : }
79 0 : }
80 :
81 0 : void storeImg(String fileName,ImageInterface<Float>& theImg)
82 : {
83 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
84 0 : LatticeExpr<Float> le(theImg);
85 0 : tmp.copyData(le);
86 0 : }
87 :
88 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
89 : const Array<Complex>& theImg)
90 : {
91 0 : PagedImage<Complex> ctmp(theImg.shape(), coord, fileName);
92 0 : ctmp.put(theImg);
93 0 : }
94 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
95 : const Array<DComplex>& theImg)
96 : {
97 0 : PagedImage<DComplex> ctmp(theImg.shape(), coord, fileName);
98 0 : ctmp.put(theImg);
99 0 : }
100 0 : void storeArrayAsImage(String fileName, const CoordinateSystem& coord,
101 : const Array<Float>& theImg)
102 : {
103 0 : PagedImage<Float> ctmp(theImg.shape(), coord, fileName);
104 0 : ctmp.put(theImg);
105 0 : }
106 : //
107 : //---------------------------------------------------------------------
108 : //
109 : // Get the time stamp for the for the current visibility integration.
110 : // Since VisTimeAverager() does not accumulate auto-correlations (it
111 : // should though!), and the VisBuffer can potentially have
112 : // auto-correlation placeholders, vb.time()(0) may not be correct (it
113 : // will be in fact zero when AC rows are present). So look for the
114 : // first timestamp of a row corresponding to an unflagged
115 : // cross-correlation.
116 : //
117 0 : Double getCurrentTimeStamp(const VisBuffer2& vb)
118 : {
119 : //LogIO os(LogOrigin("Utils", "getCurrentTimeStamp", WHERE));
120 :
121 0 : Int N=vb.nRows();
122 0 : for(Int i=0;i<N;i++)
123 : {
124 0 : if ((!vb.flagRow()(i)) && (vb.antenna1()(i) != vb.antenna2()(i)))
125 0 : return vb.time()(i);
126 : }
127 : //os << "No unflagged row found! This is a major problem/bug" << LogIO::WARN;
128 0 : return 0.0;
129 : }
130 : //
131 : //---------------------------------------------------------------------
132 : // Compute the Parallactic Angle for the give VisBuffer
133 : //
134 0 : void getHADec(MeasurementSet& ms, const VisBuffer2& vb,
135 : Double& HA, Double& RA, Double& Dec)
136 : {
137 0 : MEpoch last;
138 0 : Double time = getCurrentTimeStamp(vb);
139 :
140 0 : Unit Second("s"), Day("d");
141 0 : MEpoch epoch(Quantity(time,Second),MEpoch::TAI);
142 0 : MPosition pos;
143 0 : MSObservationColumns msoc(ms.observation());
144 0 : String ObsName=msoc.telescopeName()(vb.arrayId()(0));
145 :
146 0 : if (!MeasTable::Observatory(pos,ObsName))
147 0 : throw(AipsError("Observatory position for "+ ObsName + " not found"));
148 :
149 0 : MeasFrame frame(epoch,pos);
150 0 : MEpoch::Convert toLAST = MEpoch::Convert(MEpoch::Ref(MEpoch::TAI,frame),
151 0 : MEpoch::Ref(MEpoch::LAST,frame));
152 0 : RA=vb.direction1()(0).getAngle().getValue()(0);
153 0 : if (RA < 0.0) RA += M_PI*2.0;
154 0 : Dec=vb.direction1()(0).getAngle().getValue()(1);
155 :
156 0 : last = toLAST(epoch);
157 0 : Double LST = last.get(Day).getValue();
158 0 : LST -= floor(LST); // Extract the fractional day
159 0 : LST *= 2*C::pi;// Convert to Raidan
160 :
161 0 : cout << "LST = " << LST << " " << LST/(2*C::pi) << endl;
162 :
163 0 : HA = LST - RA;
164 0 : }
165 : //
166 : //---------------------------------------------------------------------
167 : // Compute the Parallactic Angle for the give VisBuffer
168 : //
169 0 : Double getPA(const vi::VisBuffer2& vb)
170 : {
171 0 : return (Double)(vb.feedPa(getCurrentTimeStamp(vb))(0));
172 : // Double pa=0;
173 : // Int n=0;
174 : // Vector<Float> antPA = vb.feed_pa(getCurrentTimeStamp(vb));
175 : // for (uInt i=0;i<antPA.nelements();i++)
176 : // {
177 : // if (!vb.msColumns().antenna().flagRow()(i))
178 : // {pa += antPA(i); n++;break;}
179 : // }
180 : // // pa = sum(antPA)/(antPA.nelements()-1);
181 : // if (n==0)
182 : // throw(AipsError("No unflagged antenna found in getPA()."));
183 : // return pa/n;
184 : }
185 : //
186 : //---------------------------------------------------------------------
187 : //
188 : // Make stokes axis, given the polarization types.
189 : //
190 0 : void makeStokesAxis(Int npol_p, Vector<String>& polType, Vector<Int>& whichStokes)
191 : {
192 : // Vector<String> polType=msc.feed().polarizationType()(0);
193 : StokesImageUtil::PolRep polRep_p;
194 0 : LogIO os(LogOrigin("Utils", "makeStokesAxis", WHERE));
195 :
196 0 : if (polType(0)!="X" && polType(0)!="Y" &&
197 0 : polType(0)!="R" && polType(0)!="L")
198 : {
199 : os << "Warning: Unknown stokes types in feed table: ["
200 0 : << polType(0) << ", " << polType(1) << "]" << endl
201 0 : << "Results open to question!" << LogIO::POST;
202 : }
203 :
204 0 : if (polType(0)=="X" || polType(0)=="Y")
205 : {
206 0 : polRep_p=StokesImageUtil::LINEAR;
207 0 : os << "Preferred polarization representation is linear" << LogIO::POST;
208 : }
209 : else
210 : {
211 0 : polRep_p=StokesImageUtil::CIRCULAR;
212 0 : os << "Preferred polarization representation is circular" << LogIO::POST;
213 : }
214 :
215 : // Vector<Int> whichStokes(npol_p);
216 0 : switch(npol_p)
217 : {
218 0 : case 1:
219 0 : whichStokes.resize(1);
220 0 : whichStokes(0)=Stokes::I;
221 0 : os << "Image polarization = Stokes I" << LogIO::POST;
222 0 : break;
223 0 : case 2:
224 0 : whichStokes.resize(2);
225 0 : whichStokes(0)=Stokes::I;
226 0 : if (polRep_p==StokesImageUtil::LINEAR)
227 : {
228 0 : whichStokes(1)=Stokes::Q;
229 0 : os << "Image polarization = Stokes I,Q" << LogIO::POST;
230 : }
231 : else
232 : {
233 0 : whichStokes(1)=Stokes::V;
234 0 : os << "Image polarization = Stokes I,V" << LogIO::POST;
235 : }
236 0 : break;
237 0 : case 3:
238 0 : whichStokes.resize(3);
239 0 : whichStokes(0)=Stokes::I;
240 0 : whichStokes(1)=Stokes::Q;
241 0 : whichStokes(1)=Stokes::U;
242 0 : os << "Image polarization = Stokes I,Q,U" << LogIO::POST;
243 0 : break;
244 0 : case 4:
245 0 : whichStokes.resize(4);
246 0 : whichStokes(0)=Stokes::I;
247 0 : whichStokes(1)=Stokes::Q;
248 0 : whichStokes(2)=Stokes::U;
249 0 : whichStokes(3)=Stokes::V;
250 0 : os << "Image polarization = Stokes I,Q,U,V" << LogIO::POST;
251 0 : break;
252 0 : default:
253 : os << LogIO::SEVERE << "Illegal number of Stokes parameters: " << npol_p
254 0 : << LogIO::POST;
255 : };
256 0 : }
257 : //
258 : //--------------------------------------------------------------------------------------------
259 : //
260 0 : Bool isVBNaN(const VisBuffer2 &vb,String& mesg)
261 : {
262 0 : IPosition ndx(3,0);
263 0 : Int N = vb.nRows();
264 0 : for(ndx(2)=0;ndx(2)<N;ndx(2)++)
265 0 : if (std::isnan(vb.visCubeModel()(ndx).real()) ||
266 0 : std::isnan(vb.visCubeModel()(ndx).imag())
267 : )
268 : {
269 0 : ostringstream os;
270 0 : os << ndx(2) << " " << vb.antenna1()(ndx(2)) << "-" << vb.antenna2()(ndx(2)) << " "
271 0 : << vb.flagCube()(ndx) << " "
272 : //<< vb.flag()(0,ndx(2)) << " "
273 0 : << vb.weight();
274 0 : mesg = os.str().c_str();
275 0 : return true;
276 : }
277 0 : return false;
278 : }
279 : //
280 : //--------------------------------------------------------------------------------------------
281 : //
282 : /////////////////////////////////////////////////////////////////////////////
283 : //
284 : // IChangeDetector - an interface class to detect changes in the VisBuffer
285 : // Exact meaning of the "change" is defined in the derived classes
286 : // (e.g. a change in the parallactic angle)
287 :
288 : // return true if a change occurs somewhere in the buffer
289 : using namespace vi;
290 0 : Bool IChangeDetector::changed(const VisBuffer2 &vb) const
291 : {
292 0 : for (rownr_t i=0;i<vb.nRows();++i)
293 0 : if (changed(vb,i)) return true;
294 0 : return false;
295 : }
296 :
297 : // return true if a change occurs somewhere in the buffer starting from row1
298 : // up to row2 (row2=-1 means up to the end of the buffer). The row number,
299 : // where the change occurs is returned in the row2 parameter
300 0 : Bool IChangeDetector::changedBuffer(const VisBuffer2 &vb, Int row1,
301 : Int &row2) const
302 : {
303 0 : if (row1<0) row1=0;
304 0 : int jrow = row2;
305 0 : if (jrow < 0) jrow = vb.nRows()-1;
306 0 : DebugAssert(jrow<(int)vb.nRows(),AipsError);
307 :
308 : // It is not important now to have a separate function for a "block"
309 : // operation. Because an appropriate caching is implemented inside
310 : // VisBuffer, changed(vb,row) can be called in the
311 : // loop without a perfomance penalty. This method returns the
312 : // first row where the change occured rather than the last unchanged
313 : // row as it was for BeamSkyJones::changedBuffer
314 :
315 0 : for (int ii=row1;ii<=jrow;++ii)
316 0 : if (changed(vb,ii)) {
317 0 : row2 = ii;
318 0 : return true;
319 : }
320 0 : return false;
321 : }
322 :
323 : // a virtual destructor to make the compiler happy
324 0 : IChangeDetector::~IChangeDetector() {}
325 :
326 : //
327 : /////////////////////////////////////////////////////////////////////////////
328 :
329 : /////////////////////////////////////////////////////////////////////////////
330 : //
331 : // ParAngleChangeDetector - a class to detect a change in the parallatic
332 : // angle
333 : //
334 :
335 : // set up the tolerance, which determines how much the position angle should
336 : // change to report the change by this class
337 0 : ParAngleChangeDetector::ParAngleChangeDetector(const Quantity &pa_tolerance)
338 0 : : pa_tolerance_p(pa_tolerance.getValue("rad")),
339 0 : last_pa_p(1000.) {} // 1000 is >> 2pi, so it is changed
340 : // after construction
341 :
342 : // Set the value of the PA tolerance
343 0 : void ParAngleChangeDetector::setTolerance(const Quantity &pa_tolerance)
344 : {
345 0 : pa_tolerance_p = (pa_tolerance.getValue("rad"));
346 0 : }
347 : // reset to the state which exist just after construction
348 0 : void ParAngleChangeDetector::reset()
349 : {
350 0 : last_pa_p=1000.; // it is >> 2pi, which would force a changed state
351 0 : }
352 :
353 : // return parallactic angle tolerance
354 0 : Quantity ParAngleChangeDetector::getParAngleTolerance() const
355 : {
356 0 : return Quantity(pa_tolerance_p,"rad");
357 : }
358 :
359 : // return true if a change occurs in the given row since the last call
360 : // of update
361 0 : Bool ParAngleChangeDetector::changed(const VisBuffer2 &vb, Int row) const
362 : {
363 0 : if (row<0) row=0;
364 : // const Double feed1_pa=vb.feed1_pa()[row];
365 0 : Double feed1_pa=getPA(vb);
366 0 : if (abs(feed1_pa-last_pa_p) > pa_tolerance_p)
367 : {
368 : // cout << "Utils: " << feed1_pa*57.295 << " " << last_pa_p*57.295 << " " << abs(feed1_pa-last_pa_p)*57.295 << " " << ttt*57.295 << " " << vb.time()(0)-4.51738e+09 << endl;
369 0 : return true;
370 : }
371 : // const Double feed2_pa=vb.feed2_pa()[row];
372 0 : Double feed2_pa = getPA(vb);
373 0 : if (abs(feed2_pa-last_pa_p) > pa_tolerance_p)
374 : {
375 : // cout << "Utils: " << feed2_pa*57.295 << " "
376 : // << last_pa_p*57.295 << " "
377 : // << abs(feed2_pa-last_pa_p)*57.295 << " " << ttt*57.295 << vb.time()(0)-4.51738e+09 <<endl;
378 0 : return true;
379 : }
380 0 : return false;
381 : }
382 :
383 : // start looking for a change from the given row of the VisBuffer
384 0 : void ParAngleChangeDetector::update(const VisBuffer2 &vb, Int row)
385 : {
386 0 : if (row<0) row=0;
387 0 : const Double feed1_pa=vb.feedPa1()[row];
388 0 : const Double feed2_pa=vb.feedPa2()[row];
389 0 : if (abs(feed1_pa-feed2_pa)>pa_tolerance_p) {
390 0 : LogIO os;
391 0 : os<<LogIO::WARN << LogOrigin("ParAngleChangeDetector","update")
392 : <<"The parallactic angle is different at different stations"
393 0 : <<LogIO::POST<<LogIO::WARN <<LogOrigin("ParAngleChangeDetector","update")
394 0 : <<"The result may be incorrect. Continuing anyway."<<LogIO::POST;
395 : }
396 0 : last_pa_p=(feed1_pa+feed2_pa)/2.;
397 0 : }
398 :
399 0 : Int getPhaseCenter(MeasurementSet& ms, MDirection& dir0, Int whichField)
400 : {
401 0 : MSFieldColumns msfc(ms.field());
402 0 : if (whichField < 0)
403 : {
404 0 : MSRange msr(ms);
405 : //
406 : // An array of shape [2,1,1]!
407 : //
408 0 : Vector<Int> fieldId;
409 0 : fieldId = msr.range(MSS::FIELD_ID).asArrayInt(RecordFieldId(0));
410 :
411 0 : Array<Double> phaseDir = msfc.phaseDir().getColumn();
412 :
413 0 : if (fieldId.nelements() <= 0)
414 0 : throw(AipsError("getPhaseCenter: No fields found in the selected MS."));
415 :
416 0 : IPosition ndx0(3,0,0,0),ndx1(3,1,0,0);
417 :
418 : Double maxRA, maxDec, RA,Dec,RA0,Dec0, dist0;
419 0 : maxRA = maxDec=0;
420 0 : for(uInt i=0;i<fieldId.nelements();i++)
421 : {
422 0 : RA = phaseDir(IPosition(3,0,0,fieldId(i)));
423 0 : Dec = phaseDir(IPosition(3,1,0,fieldId(i)));
424 0 : maxRA += RA; maxDec += Dec;
425 : }
426 0 : RA0=maxRA/fieldId.nelements(); Dec0=maxDec/fieldId.nelements();
427 :
428 0 : dist0=10000;
429 0 : Int field0=0;
430 0 : for(uInt i=0;i<fieldId.nelements();i++)
431 : {
432 0 : RA = RA0-phaseDir(IPosition(3,0,0,fieldId(i)));
433 0 : Dec = Dec0-phaseDir(IPosition(3,1,0,fieldId(i)));
434 0 : Double dist=sqrt(RA*RA + Dec*Dec);
435 0 : if (dist < dist0) {field0=fieldId(i);dist0=dist;};
436 : }
437 0 : dir0=msfc.phaseDirMeasCol()(field0)(IPosition(1,0));
438 : //dir0=msfc.phaseDirMeasCol()(6)(IPosition(1,0));
439 0 : return field0;
440 : }
441 : else
442 : {
443 0 : dir0=msfc.phaseDirMeasCol()(whichField)(IPosition(1,0));
444 0 : return whichField;
445 : }
446 : }
447 : //
448 : //------------------------------------------------------------------
449 : //
450 0 : Bool findMaxAbsLattice(const ImageInterface<Float>& lattice,
451 : Float& maxAbs,IPosition& posMaxAbs)
452 : {
453 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
454 0 : maxAbs=0.0;
455 :
456 0 : const IPosition tileShape = lattice.niceCursorShape();
457 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
458 : {
459 0 : RO_LatticeIterator<Float> li(lattice, ls);
460 0 : for(li.reset();!li.atEnd();li++)
461 : {
462 0 : IPosition posMax=li.position();
463 0 : IPosition posMin=li.position();
464 0 : Float maxVal=0.0;
465 0 : Float minVal=0.0;
466 :
467 0 : minMax(minVal, maxVal, posMin, posMax, li.cursor());
468 0 : if((maxVal)>(maxAbs))
469 : {
470 0 : maxAbs=maxVal;
471 0 : posMaxAbs=li.position();
472 0 : posMaxAbs(0)=posMax(0);
473 : }
474 : }
475 : }
476 :
477 0 : return true;
478 : };
479 : //
480 : //------------------------------------------------------------------
481 : //
482 0 : Bool findMaxAbsLattice(const ImageInterface<Float>& masklat,
483 : const Lattice<Float>& lattice,
484 : Float& maxAbs,IPosition& posMaxAbs,
485 : Bool flip)
486 : {
487 :
488 0 : AlwaysAssert(masklat.shape()==lattice.shape(), AipsError);
489 :
490 0 : Array<Float> msk;
491 :
492 0 : posMaxAbs = IPosition(lattice.shape().nelements(), 0);
493 0 : maxAbs=0.0;
494 : //maxAbs=-1.0e+10;
495 0 : const IPosition tileShape = lattice.niceCursorShape();
496 0 : TiledLineStepper ls(lattice.shape(), tileShape, 0);
497 0 : TiledLineStepper lsm(masklat.shape(), tileShape, 0);
498 : {
499 0 : RO_LatticeIterator<Float> li(lattice, ls);
500 0 : RO_LatticeIterator<Float> lim(masklat, lsm);
501 0 : for(li.reset(),lim.reset();!li.atEnd();li++,lim++)
502 : {
503 0 : IPosition posMax=li.position();
504 0 : IPosition posMin=li.position();
505 0 : Float maxVal=0.0;
506 0 : Float minVal=0.0;
507 :
508 0 : msk = lim.cursor();
509 0 : if(flip) msk = (Float)1.0 - msk;
510 :
511 : //minMaxMasked(minVal, maxVal, posMin, posMax, li.cursor(),lim.cursor());
512 0 : minMaxMasked(minVal, maxVal, posMin, posMax, li.cursor(),msk);
513 :
514 :
515 0 : if((maxVal)>(maxAbs))
516 : {
517 0 : maxAbs=maxVal;
518 0 : posMaxAbs=li.position();
519 0 : posMaxAbs(0)=posMax(0);
520 : }
521 : }
522 : }
523 :
524 0 : return true;
525 : }
526 : //
527 : //---------------------------------------------------------------
528 : //Rotate a complex array using a the given coordinate system and the
529 : //angle in radians. Default interpolation method is "CUBIC".
530 : //Axeses corresponding to Linear coordinates in the given
531 : //CoordinateSystem object are rotated. Rotation is done using
532 : //ImageRegrid object, about the pixel given by (N+1)/2 where N is
533 : //the number of pixels along the axis.
534 : //
535 0 : void SynthesisUtils::rotateComplexArray(LogIO& logio, Array<Complex>& inArray,
536 : CoordinateSystem& inCS,
537 : Array<Complex>& outArray,
538 : Double dAngleRad,
539 : String interpMethod,
540 : Bool modifyInCS)
541 : {
542 : // logio << LogOrigin("SynthesisUtils", "rotateComplexArray")
543 : // << "Rotating CF using " << interpMethod << " interpolation."
544 : // << LogIO::POST;
545 : (void)logio;
546 : //
547 : // If no rotation required, just copy the inArray to outArray.
548 : //
549 : // if (abs(dAngleRad) < 0.1)
550 :
551 : // IPosition tt;
552 : // inCS.list(logio,MDoppler::RADIO,tt,tt);
553 :
554 0 : if (abs(dAngleRad) == 0.0)
555 : {
556 0 : outArray.reference(inArray);
557 : // outArray.assign(inArray);
558 0 : return;
559 : }
560 : //
561 : // Re-grid inImage onto outImage
562 : //
563 0 : Vector<Int> pixelAxes;
564 0 : Int coordInd = -1;
565 : // Extract LINRAR coords from inCS.
566 : // Extract axes2
567 :
568 0 : if(modifyInCS){
569 0 : Vector<Double> refPix = inCS.referencePixel();
570 0 : refPix(0) = (Int)((inArray.shape()(0))/2.0);
571 0 : refPix(1) = (Int)((inArray.shape()(1))/2.0);
572 0 : inCS.setReferencePixel(refPix);
573 : }
574 :
575 0 : coordInd = inCS.findCoordinate(Coordinate::LINEAR);
576 0 : Bool haveLinear = true;
577 :
578 0 : if(coordInd == -1){ // no linear coordinate found, look for DIRECTION instead
579 0 : coordInd = inCS.findCoordinate(Coordinate::DIRECTION);
580 0 : haveLinear = false;
581 : }
582 :
583 0 : pixelAxes=inCS.pixelAxes(coordInd);
584 0 : IPosition axes2(pixelAxes);
585 : // Set linear transformation matrix in inCS.
586 : // CoordinateSystem outCS =
587 : // ImageRegrid<Complex>::makeCoordinateSystem (logio, outCS, inCS, axes2);
588 :
589 0 : CoordinateSystem outCS(inCS);
590 :
591 0 : Matrix<Double> xf = outCS.coordinate(coordInd).linearTransform();
592 0 : Matrix<Double> rotm(2,2);
593 0 : rotm(0,0) = cos(dAngleRad); rotm(0,1) = sin(dAngleRad);
594 0 : rotm(1,0) = -rotm(0,1); rotm(1,1) = rotm(0,0);
595 :
596 : // Double s = sin(dAngleRad);
597 : // Double c = cos(dAngleRad);
598 : // rotm(0,0) = c; rotm(0,1) = s;
599 : // rotm(1,0) = -s; rotm(1,1) = c;
600 :
601 : // Create new linear transform matrix
602 0 : Matrix<Double> xform(2,2);
603 0 : xform(0,0) = rotm(0,0)*xf(0,0)+rotm(0,1)*xf(1,0);
604 0 : xform(0,1) = rotm(0,0)*xf(0,1)+rotm(0,1)*xf(1,1);
605 0 : xform(1,0) = rotm(1,0)*xf(0,0)+rotm(1,1)*xf(1,0);
606 0 : xform(1,1) = rotm(1,0)*xf(0,1)+rotm(1,1)*xf(1,1);
607 :
608 0 : if(haveLinear){
609 0 : LinearCoordinate linCoords = outCS.linearCoordinate(coordInd);
610 0 : linCoords.setLinearTransform(xform);
611 0 : outCS.replaceCoordinate(linCoords, coordInd);
612 : }
613 : else{
614 0 : DirectionCoordinate dirCoords = outCS.directionCoordinate(coordInd);
615 0 : dirCoords.setLinearTransform(xform);
616 0 : outCS.replaceCoordinate(dirCoords, coordInd);
617 : }
618 :
619 0 : outArray.resize(inArray.shape());
620 0 : outArray.set(0);
621 : //
622 : // Make an image out of inArray and inCS --> inImage
623 : //
624 : // TempImage<Complex> inImage(inArray.shape(), inCS);
625 : {
626 0 : TempImage<Float> inImage(inArray.shape(),inCS);
627 0 : TempImage<Float> outImage(outArray.shape(), outCS);
628 0 : ImageRegrid<Float> ir;
629 0 : Interpolate2D::Method interpolationMethod = Interpolate2D::stringToMethod(interpMethod);
630 : //------------------------------------------------------------------------
631 : // Rotated the real part
632 : //
633 0 : inImage.copyData(LatticeExpr<Float>(real(ArrayLattice<Complex>(inArray))));
634 0 : outImage.set(0.0);
635 :
636 0 : ir.regrid(outImage, interpolationMethod, axes2, inImage);
637 0 : setReal(outArray,outImage.get());
638 : //------------------------------------------------------------------------
639 : // Rotated the imaginary part
640 : //
641 0 : inImage.copyData(LatticeExpr<Float>(imag(ArrayLattice<Complex>(inArray))));
642 0 : outImage.set(0.0);
643 :
644 0 : ir.regrid(outImage, interpolationMethod, axes2, inImage);
645 0 : setImag(outArray,outImage.get());
646 : }
647 : }
648 : //
649 : //---------------------------------------------------------------
650 : //
651 0 : void SynthesisUtils::findLatticeMax(const ImageInterface<Complex>& lattice,
652 : Vector<Float>& maxAbs,
653 : Vector<IPosition>& posMaxAbs)
654 : {
655 0 : IPosition lshape(lattice.shape());
656 0 : IPosition ndx(lshape);
657 0 : Int nPol=lshape(2);
658 0 : posMaxAbs.resize(nPol);
659 0 : for(Int i=0;i<nPol;i++)
660 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
661 0 : maxAbs.resize(nPol);
662 0 : ndx=0;
663 :
664 0 : for(Int s2=0;s2<lshape(2);s2++)
665 0 : for(Int s3=0;s3<lshape(3);s3++)
666 : {
667 0 : ndx(2) = s2; ndx(3)=s3;
668 : {
669 : //
670 : // Locate the pixel with the peak value. That's the
671 : // origin in pixel co-ordinates.
672 : //
673 0 : maxAbs(s2)=0;
674 0 : posMaxAbs(s2) = 0;
675 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
676 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
677 0 : if (abs(lattice(ndx)) > maxAbs(s2))
678 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
679 : }
680 : }
681 0 : }
682 : //
683 : //---------------------------------------------------------------
684 : //
685 0 : void SynthesisUtils::findLatticeMax(const Array<Complex>& lattice,
686 : Vector<Float>& maxAbs,
687 : Vector<IPosition>& posMaxAbs)
688 : {
689 0 : IPosition lshape(lattice.shape());
690 0 : IPosition ndx(lshape);
691 0 : Int nPol=lshape(2);
692 0 : posMaxAbs.resize(nPol);
693 0 : for(Int i=0;i<nPol;i++)
694 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
695 0 : maxAbs.resize(nPol);
696 0 : ndx=0;
697 :
698 0 : for(Int s2=0;s2<lshape(2);s2++)
699 0 : for(Int s3=0;s3<lshape(3);s3++)
700 : {
701 0 : ndx(2) = s2; ndx(3)=s3;
702 : {
703 : //
704 : // Locate the pixel with the peak value. That's the
705 : // origin in pixel co-ordinates.
706 : //
707 0 : maxAbs(s2)=0;
708 0 : posMaxAbs(s2) = 0;
709 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
710 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
711 0 : if (abs(lattice(ndx)) > maxAbs(s2))
712 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
713 : }
714 : }
715 0 : }
716 : //
717 : //---------------------------------------------------------------
718 : //
719 0 : void SynthesisUtils::findLatticeMax(const ImageInterface<Float>& lattice,
720 : Vector<Float>& maxAbs,
721 : Vector<IPosition>& posMaxAbs)
722 : {
723 0 : IPosition lshape(lattice.shape());
724 0 : IPosition ndx(lshape);
725 0 : Int nPol=lshape(2);
726 0 : posMaxAbs.resize(nPol);
727 0 : for(Int i=0;i<nPol;i++)
728 0 : posMaxAbs(i)=IPosition(lattice.shape().nelements(), 0);
729 0 : maxAbs.resize(nPol);
730 0 : ndx=0;
731 :
732 0 : for(Int s2=0;s2<lshape(2);s2++)
733 0 : for(Int s3=0;s3<lshape(3);s3++)
734 : {
735 0 : ndx(2) = s2; ndx(3)=s3;
736 : {
737 : //
738 : // Locate the pixel with the peak value. That's the
739 : // origin in pixel co-ordinates.
740 : //
741 0 : maxAbs(s2)=0;
742 0 : posMaxAbs(s2) = 0;
743 0 : for(ndx(1)=0;ndx(1)<lshape(1);ndx(1)++)
744 0 : for(ndx(0)=0;ndx(0)<lshape(0);ndx(0)++)
745 0 : if (abs(lattice(ndx)) > maxAbs(s2))
746 0 : {posMaxAbs(s2) = ndx;maxAbs(s2)=abs(lattice(ndx));}
747 : }
748 : }
749 0 : }
750 : //
751 : //---------------------------------------------------------------
752 : // Get the value of the named variable from ~/.aipsrc (or ~/.casarc)
753 : // or from a env. variable (in this precidence order).
754 : //
755 : template <class T>
756 0 : T SynthesisUtils::getenv(const char *name, const T defaultVal)
757 : {
758 0 : T val=defaultVal;
759 0 : stringstream defaultStr;
760 0 : defaultStr << defaultVal;
761 : {
762 0 : char *valStr=NULL;
763 0 : std::string tt(name);
764 : unsigned long pos;
765 0 : while((pos=tt.find(".")) != tt.npos)
766 0 : tt.replace(pos, 1, "_");
767 :
768 0 : if ((valStr = std::getenv(tt.c_str())) != NULL)
769 : {
770 0 : stringstream toT2(valStr);
771 0 : toT2 >> val;
772 : }
773 : }
774 : // If environment variable was not found (val remained set to the
775 : // defaulVal), look in ~/.aipsrc.
776 0 : if (val==defaultVal)
777 : {
778 0 : uint handle = Aipsrc::registerRC(name, defaultStr.str().c_str());
779 0 : String strVal = Aipsrc::get(handle);
780 0 : stringstream toT(strVal);
781 0 : toT >> val;
782 : }
783 0 : return val;
784 : }
785 : template
786 : Int SynthesisUtils::getenv(const char *name, const Int defaultVal);
787 : template
788 : Bool SynthesisUtils::getenv(const char *name, const Bool defaultVal);
789 : template
790 : Float SynthesisUtils::getenv(const char *name, const Float defaultVal);
791 : template
792 : double SynthesisUtils::getenv(const char *name, const double defaultVal);
793 : template
794 : String SynthesisUtils::getenv(const char *name, const String defaultVal);
795 :
796 0 : Float SynthesisUtils::libreSpheroidal(Float nu)
797 : {
798 : Double top, bot, nuend, delnusq;
799 : uInt part;
800 0 : if (nu <= 0) return 1.0;
801 : else
802 0 : if (nu >= 1.0)
803 0 : return 0.0;
804 : else
805 : {
806 0 : uInt np = 5;
807 0 : uInt nq = 3;
808 0 : Matrix<Double> p(np, 2);
809 0 : Matrix<Double> q(nq, 2);
810 0 : p(0,0) = 8.203343e-2;
811 0 : p(1,0) = -3.644705e-1;
812 0 : p(2,0) = 6.278660e-1;
813 0 : p(3,0) = -5.335581e-1;
814 0 : p(4,0) = 2.312756e-1;
815 0 : p(0,1) = 4.028559e-3;
816 0 : p(1,1) = -3.697768e-2;
817 0 : p(2,1) = 1.021332e-1;
818 0 : p(3,1) = -1.201436e-1;
819 0 : p(4,1) = 6.412774e-2;
820 0 : q(0,0) = 1.0000000e0;
821 0 : q(1,0) = 8.212018e-1;
822 0 : q(2,0) = 2.078043e-1;
823 0 : q(0,1) = 1.0000000e0;
824 0 : q(1,1) = 9.599102e-1;
825 0 : q(2,1) = 2.918724e-1;
826 :
827 0 : part = 0;
828 0 : nuend = 0.0;
829 0 : if ((nu >= 0.0) && (nu < 0.75))
830 : {
831 0 : part = 0;
832 0 : nuend = 0.75;
833 : }
834 0 : else if ((nu >= 0.75) && (nu <= 1.00))
835 : {
836 0 : part = 1;
837 0 : nuend = 1.0;
838 : }
839 :
840 0 : top = p(0,part);
841 0 : delnusq = pow(nu,2.0) - pow(nuend,2.0);
842 0 : for (uInt k=1; k<np; k++)
843 0 : top += p(k, part) * pow(delnusq, (Double)k);
844 :
845 0 : bot = q(0, part);
846 0 : for (uInt k=1; k<nq; k++)
847 0 : bot += q(k,part) * pow(delnusq, (Double)k);
848 :
849 0 : if (bot != 0.0) return (top/bot);
850 0 : else return 0.0;
851 : }
852 : }
853 0 : Double SynthesisUtils::getRefFreq(const VisBuffer2& /*vb*/)
854 : {
855 0 : throw(AipsError("SynthesisUtils::getRefFreq() depricated due to VI2/VB2 move"));
856 : // return max((vb.getVi()->ms())//vb.msColumns()
857 : // .spectralWindow().chanFreq().getColumn());
858 : }
859 :
860 0 : void SynthesisUtils::makeFTCoordSys(const CoordinateSystem& coords,
861 : const Int& convSize,
862 : const Vector<Double>& ftRef,
863 : CoordinateSystem& ftCoords)
864 : {
865 : Int directionIndex;
866 :
867 0 : ftCoords = coords;
868 0 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
869 : // The following line follows the (lame) logic that if a
870 : // DIRECTION axis was not found, the coordinate system must be of
871 : // the FT domain already
872 0 : if (directionIndex == -1) return;
873 :
874 0 : DirectionCoordinate dc;//=coords.directionCoordinate(directionIndex);
875 : // AlwaysAssert(directionIndex>=0, AipsError);
876 0 : dc=coords.directionCoordinate(directionIndex);
877 0 : Vector<Bool> axes(2); axes(0)=axes(1)=true;//axes(2)=true;
878 0 : Vector<Int> shape(2,convSize);
879 :
880 0 : Vector<Double>ref(4);
881 0 : ref(0)=ref(1)=ref(2)=ref(3)=0;
882 0 : dc.setReferencePixel(ref);
883 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
884 0 : Vector<Double> refVal;
885 0 : refVal=ftdc->referenceValue();
886 : // refVal(0)=refVal(1)=0;
887 : // ftdc->setReferenceValue(refVal);
888 0 : ref(0)=ftRef(0);
889 0 : ref(1)=ftRef(1);
890 0 : ftdc->setReferencePixel(ref);
891 :
892 0 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
893 : // LogIO logio;
894 : // IPosition tt;
895 : // coords.list(logio,MDoppler::RADIO,tt,tt);
896 : // ftCoords.list(logio,MDoppler::RADIO,tt,tt);
897 :
898 0 : delete ftdc; ftdc=0;
899 : }
900 :
901 : //
902 : // Given a list of Spw,MinFreq,MaxFreq,FreqStep (the output product
903 : // of MSSelection), expand the list to a list of channel freqs. and
904 : // conjugate freqs. per SPW.
905 : //
906 0 : void SynthesisUtils::expandFreqSelection(const Matrix<Double>& freqSelection,
907 : Matrix<Double>& expandedFreqList,
908 : Matrix<Double>& expandedConjFreqList)
909 : {
910 0 : Int nSpw = freqSelection.shape()(0), maxSlots=0;
911 : Double freq;
912 :
913 0 : for (Int s=0;s<nSpw;s++)
914 0 : maxSlots=max(maxSlots,SynthesisUtils::nint((freqSelection(s,2)-freqSelection(s,1))/freqSelection(s,3))+1);
915 :
916 0 : expandedFreqList.resize(nSpw,maxSlots);
917 0 : expandedConjFreqList.resize(nSpw,maxSlots);
918 :
919 0 : for (Int s=0,cs=(nSpw-1);s<nSpw;s++,cs--)
920 0 : for (Int i=0,ci=(maxSlots-1);i<maxSlots;i++,ci--)
921 : {
922 0 : freq = freqSelection(s,1)+i*freqSelection(s,3);
923 0 : expandedFreqList(s,i) = (freq <= freqSelection(s,2)) ? freq : 0;
924 0 : freq = freqSelection(cs,2) - ci*freqSelection(cs,3);
925 0 : expandedConjFreqList(s,ci) = (freq >= freqSelection(cs,1)) ? freq : 0;
926 : }
927 0 : }
928 :
929 : //
930 : // The result will be in-place in c1
931 : //
932 : template
933 : void SynthesisUtils::libreConvolver(Array<Complex>& c1, const Array<Complex>& c2);
934 :
935 :
936 : template <class T>
937 0 : void SynthesisUtils::libreConvolver(Array<T>& c1, const Array<T>& c2)
938 : {
939 0 : Array<T> c2tmp;
940 0 : c2tmp.assign(c2);
941 :
942 0 : if (c1.shape().product() > c2tmp.shape().product())
943 0 : c2tmp.resize(c1.shape(),true);
944 : else
945 0 : c1.resize(c2tmp.shape(),true);
946 :
947 :
948 0 : ArrayLattice<T> c2tmp_lat(c2tmp), c1_lat(c1);
949 :
950 0 : LatticeFFT::cfft2d(c1_lat,false);
951 0 : LatticeFFT::cfft2d(c2tmp_lat,false);
952 : //cerr << "########## " << c1.shape() << " " << c2tmp.shape() << endl;
953 0 : c1 = sqrt(c1);
954 0 : c2tmp=sqrt(c2tmp);
955 0 : c1 *= conj(c2tmp);
956 0 : LatticeFFT::cfft2d(c1_lat);
957 0 : }
958 :
959 :
960 0 : Double SynthesisUtils::nearestValue(const Vector<Double>& list, const Double& val, Int& index)
961 : {
962 0 : Vector<Double> diff = fabs(list - val);
963 0 : Double minVal=1e20;
964 0 : Int i=0;
965 :
966 0 : for (index=0;index<(Int)diff.nelements();index++)
967 0 : if (diff[index] < minVal) {minVal=diff[index];i=index;}
968 0 : index=i;
969 0 : return list(index);
970 :
971 : // The algorithm below has a N*log(N) cost.
972 : //
973 : // Bool dummy;
974 : // Sort sorter(diff.getStorage(dummy), sizeof(Double));
975 : // sorter.sortKey((uInt)0,TpDouble);
976 : // Int nch=list.nelements();
977 : // Vector<uInt> sortIndx;
978 : // sorter.sort(sortIndx, nch);
979 :
980 : // index=sortIndx(0);
981 : // return list(index);
982 :
983 :
984 : // Works only for regularly samplied list
985 : //
986 : // Int ndx=min(freqValues_p.nelements()-1,max(0,SynthesisUtils::nint((freqVal-freqValues_p[0])/freqValIncr_p)));
987 : // return ndx;
988 : }
989 :
990 : template <class T>
991 0 : T SynthesisUtils::stdNearestValue(const vector<T>& list, const T& val, Int& index)
992 : {
993 : // auto const it = std::lower_bound(list.begin(), list.end(), val);
994 : // if (it == list.begin()) return list[0];
995 : // else return list[*(it-1)];
996 :
997 0 : vector<T> diff=list;
998 0 : for (uInt i=0;i<list.size();i++)
999 0 : diff[i] = fabs(list[i] - val);
1000 :
1001 0 : T minVal=std::numeric_limits<T>::max();//1e20;
1002 0 : Int i=0;
1003 :
1004 0 : for (index=0;index<(Int)diff.size();index++)
1005 0 : if (diff[index] < minVal) {minVal=diff[index];i=index;}
1006 0 : index=i;
1007 0 : return list[index];
1008 : }
1009 :
1010 0 : CoordinateSystem SynthesisUtils::makeUVCoords(CoordinateSystem& imageCoordSys,
1011 : IPosition& shape)
1012 : {
1013 0 : CoordinateSystem FTCoords = imageCoordSys;
1014 :
1015 0 : Int dirIndex=FTCoords.findCoordinate(Coordinate::DIRECTION);
1016 0 : DirectionCoordinate dc=imageCoordSys.directionCoordinate(dirIndex);
1017 0 : Vector<Bool> axes(2); axes=true;
1018 0 : Vector<Int> dirShape(2); dirShape(0)=shape(0);dirShape(1)=shape(1);
1019 0 : Coordinate* FTdc=dc.makeFourierCoordinate(axes,dirShape);
1020 0 : FTCoords.replaceCoordinate(*FTdc,dirIndex);
1021 0 : delete FTdc;
1022 :
1023 0 : return FTCoords;
1024 : }
1025 :
1026 0 : Vector<Int> SynthesisUtils::mapSpwIDToDDID(const VisBuffer2& vb, const Int& spwID)
1027 : {
1028 0 : Vector<Int> ddidList;
1029 : //Int nDDRows = vb.msColumns().dataDescription().nrow();
1030 0 : Int nDDRows = (vb.ms()).dataDescription().nrow();
1031 0 : for (Int i=0; i<nDDRows; i++)
1032 : {
1033 0 : if((vb.subtableColumns()).dataDescription().spectralWindowId().get(i) == spwID)
1034 : {
1035 0 : Int n=ddidList.nelements();
1036 0 : ddidList.resize(n+1,true);
1037 0 : ddidList(n) = i;
1038 : }
1039 : }
1040 0 : return ddidList;
1041 : }
1042 :
1043 0 : Vector<Int> SynthesisUtils::mapSpwIDToPolID(const VisBuffer2& vb, const Int& spwID)
1044 : {
1045 0 : Vector<Int> polIDList, ddIDList;
1046 0 : ddIDList = SynthesisUtils::mapSpwIDToDDID(vb, spwID);
1047 0 : Int n=ddIDList.nelements();
1048 0 : polIDList.resize(n);
1049 0 : for (Int i=0; i<n; i++)
1050 0 : polIDList(i) = (vb.subtableColumns()).dataDescription().polarizationId().get(ddIDList(i));
1051 :
1052 0 : return polIDList;
1053 : }
1054 :
1055 :
1056 : //
1057 : // Calcualte the BLC, TRC of the intersection of two rectangles (courtesy U.Rau)
1058 : //
1059 0 : void SynthesisUtils::calcIntersection(const Int blc1[2], const Int trc1[2],
1060 : const Float blc2[2], const Float trc2[2],
1061 : Float blc[2], Float trc[2])
1062 : {
1063 : // cerr << blc1[0] << " " << blc1[1] << " " << trc1[0] << " " << trc1[1] << " " << blc2[0] << " " << blc2[1] << " " << trc2[0] << " " << trc2[1] << endl;
1064 : Float dblc, dtrc;
1065 0 : for (Int i=0;i<2;i++)
1066 : {
1067 0 : dblc = blc2[i] - blc1[i];
1068 0 : dtrc = trc2[i] - trc1[i];
1069 :
1070 0 : if ((dblc >= 0) and (dtrc >= 0))
1071 : {
1072 0 : blc[i] = blc1[i] + dblc;
1073 0 : trc[i] = trc2[i] - dtrc;
1074 : }
1075 0 : else if ((dblc >= 0) and (dtrc < 0))
1076 : {
1077 0 : blc[i] = blc1[i] + dblc;
1078 0 : trc[i] = trc1[i] + dtrc;
1079 : }
1080 0 : else if ((dblc < 0) and (dtrc >= 0))
1081 : {
1082 0 : blc[i] = blc2[i] - dblc;
1083 0 : trc[i] = trc2[i] - dtrc;
1084 : }
1085 : else
1086 : {
1087 0 : blc[i] = blc2[i] - dblc;
1088 0 : trc[i] = trc1[i] + dtrc;
1089 : }
1090 : }
1091 0 : }
1092 : //
1093 : // Check if the two rectangles interset (courtesy U.Rau).
1094 : //
1095 0 : Bool SynthesisUtils::checkIntersection(const Int blc1[2], const Int trc1[2], const Float blc2[2], const Float trc2[2])
1096 : {
1097 : // blc1[2] = {xmin1, ymin1};
1098 : // blc2[2] = {xmin2, ymin2};
1099 : // trc1[2] = {xmax1, ymax1};
1100 : // trc2[2] = {xmax2, ymax2};
1101 :
1102 0 : if ((blc1[0] > trc2[0]) || (trc1[0] < blc2[0]) || (blc1[1] > trc2[1]) || (trc1[1] < blc2[1]))
1103 0 : return false;
1104 : else
1105 0 : return true;
1106 : // def checkintersect( xmin1, ymin1, xmax1, ymax1, xmin2, ymin2, xmax2, ymax2 ):
1107 : // if xmin1 > xmax2 or xmax1 < xmin2 or ymin1 > ymax2 or ymax1 < ymin2 :
1108 : // return false
1109 : // else :
1110 : // return true
1111 : }
1112 :
1113 : // template<class Iterator>
1114 : // Iterator SynthesisUtils::Unique(Iterator first, Iterator last)
1115 : // {
1116 : // while (first != last)
1117 : // {
1118 : // Iterator next(first);
1119 : // last = std::remove(++next, last, *first);
1120 : // first = next;
1121 : // }
1122 :
1123 : // return last;
1124 : // }
1125 :
1126 0 : String SynthesisUtils::mjdToString(casacore::Time& time)
1127 : {
1128 0 : String tStr;
1129 0 : tStr = String::toString(time.year()) + "/" +
1130 0 : String::toString(time.month()) + "/" +
1131 0 : String::toString(time.dayOfMonth()) + "/" +
1132 0 : String::toString(time.hours()) + ":" +
1133 0 : String::toString(time.minutes()) + ":";
1134 0 : ostringstream fsec;
1135 0 : fsec << setprecision(2) << time.dseconds();
1136 0 : tStr = tStr + String(fsec.str());
1137 : // String::toString(time.dseconds());
1138 0 : return tStr;
1139 : }
1140 :
1141 : template <class Iterator>
1142 0 : Iterator SynthesisUtils::Unique(Iterator first, Iterator last)
1143 : {
1144 0 : while (first != last)
1145 : {
1146 0 : Iterator next(first);
1147 0 : last = std::remove(++next, last, *first);
1148 0 : first = next;
1149 : }
1150 :
1151 0 : return last;
1152 : }
1153 :
1154 0 : void SynthesisUtils::showCS(const CoordinateSystem& cs, std::ostream &os, const String& msg)
1155 : {
1156 0 : LogIO log_l(LogOrigin("SynthesisUtils","showCS"));
1157 0 : IPosition dummy;
1158 0 : Vector<String> csList;
1159 0 : if (msg!="")
1160 0 : os << "CoordSys: ";
1161 0 : csList = cs.list(log_l,MDoppler::RADIO,dummy,dummy);
1162 0 : os << csList << std::endl;
1163 0 : }
1164 0 : const casacore::Array<Complex> SynthesisUtils::getCFPixels(const casacore::String& Dir,
1165 : const casacore::String& fileName)
1166 : {
1167 : try
1168 : {
1169 0 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1170 0 : return thisCF.get();
1171 : }
1172 0 : catch (AipsError &x)
1173 : {
1174 0 : LogIO log_l(LogOrigin("SynthesisUtils","getCFPixels"));
1175 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1176 : }
1177 0 : return casacore::Array<Complex>(); // Just to keep the complier happy. Program control should never get here.
1178 : }
1179 :
1180 0 : void SynthesisUtils::putCFPixels(const casacore::String& Dir,
1181 : const casacore::String& fileName,
1182 : const casacore::Array<Complex>& srcpix)
1183 : {
1184 : try
1185 : {
1186 0 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1187 0 : return thisCF.put(srcpix);
1188 : }
1189 0 : catch (AipsError &x)
1190 : {
1191 0 : LogIO log_l(LogOrigin("SynthesisUtils","putCFPixels"));
1192 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1193 : }
1194 : }
1195 :
1196 0 : const casacore::IPosition SynthesisUtils::getCFShape(const casacore::String& Dir,
1197 : const casacore::String& fileName)
1198 : {
1199 : try
1200 : {
1201 0 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1202 0 : return thisCF.shape();
1203 : }
1204 0 : catch (AipsError &x)
1205 : {
1206 0 : LogIO log_l(LogOrigin("SynthesisUtils","getCFShape"));
1207 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1208 : }
1209 0 : return casacore::IPosition(); // Just to keep the complier happy. Program control should never get here.
1210 : }
1211 :
1212 0 : casacore::TableRecord SynthesisUtils::getCFParams(const casacore::String& Dir,
1213 : const casacore::String& fileName,
1214 : casacore::IPosition& cfShape,
1215 : casacore::Array<Complex>& pixelBuffer,
1216 : casacore::CoordinateSystem& coordSys,
1217 : casacore::Double& sampling,
1218 : casacore::Double& paVal,
1219 : casacore::Int& xSupport, casacore::Int& ySupport,
1220 : casacore::Double& fVal, casacore::Double& wVal, casacore::Int& mVal,
1221 : casacore::Double& conjFreq, casacore::Int& conjPoln,
1222 : casacore::Bool loadPixels,
1223 : casacore::Bool loadMiscInfo)
1224 : {
1225 : try
1226 : {
1227 : //casacore::Table tThisCF(Dir+'/'+fileName);
1228 0 : casacore::PagedImage<casacore::Complex> thisCF(Dir+'/'+fileName);
1229 0 : cfShape = thisCF.shape();
1230 0 : if (loadPixels) pixelBuffer.assign(thisCF.get());
1231 0 : casacore::TableRecord miscinfo;
1232 0 : if (loadMiscInfo)
1233 : {
1234 0 : miscinfo= thisCF.miscInfo();
1235 :
1236 0 : miscinfo.get("ParallacticAngle", paVal);
1237 0 : miscinfo.get("MuellerElement", mVal);
1238 0 : miscinfo.get("WValue", wVal);
1239 0 : miscinfo.get("Xsupport", xSupport);
1240 0 : miscinfo.get("Ysupport", ySupport);
1241 0 : miscinfo.get("Sampling", sampling);
1242 0 : miscinfo.get("ConjFreq", conjFreq);
1243 0 : miscinfo.get("ConjPoln", conjPoln);
1244 0 : casacore::Int index= thisCF.coordinates().findCoordinate(casacore::Coordinate::SPECTRAL);
1245 0 : coordSys = thisCF.coordinates();
1246 0 : casacore::SpectralCoordinate spCS = coordSys.spectralCoordinate(index);
1247 0 : fVal=static_cast<Float>(spCS.referenceValue()(0));
1248 : }
1249 0 : return miscinfo;
1250 : }
1251 0 : catch(AipsError& x)
1252 : {
1253 0 : throw(AipsError(String("Error in SynthesisUtils::getCFParams(): ")
1254 0 : +x.getMesg()));
1255 : }
1256 : };
1257 :
1258 :
1259 0 : void SynthesisUtils::rotate2(const double& actualPA, CFCell& baseCFC, CFCell& cfc, const double& rotAngleIncr)
1260 : {
1261 0 : LogIO log_l(LogOrigin("SynthesisUtils", "rotate2"));
1262 :
1263 : // // If the A-Term is a No-Op, the resulting CF is rotationally
1264 : // // symmetric.
1265 : // if (isNoOp()) return;
1266 :
1267 : (void)baseCFC;
1268 :
1269 : //double actualPA = getPA(vb);
1270 0 : double currentCFPA = cfc.pa_p.getValue("rad");
1271 : //Double baseCFCPA=baseCFC.pa_p.getValue("rad");
1272 :
1273 0 : double dPA = currentCFPA-actualPA;
1274 :
1275 0 : if (fabs(dPA) > fabs(rotAngleIncr))
1276 : {
1277 0 : casacore::Array<TT> inData;
1278 : //inData.assign(*baseCFC.getStorage());
1279 : //dPA = baseCFCPA-actualPA;
1280 0 : dPA = currentCFPA-actualPA;
1281 0 : inData.assign(*cfc.getStorage());
1282 : try
1283 : {
1284 0 : SynthesisUtils::rotateComplexArray(log_l, inData, cfc.coordSys_p,
1285 0 : *cfc.getStorage(),
1286 : dPA);//,"LINEAR");
1287 : // currentCFPA-actualPA);//,"LINEAR");
1288 : }
1289 0 : catch (AipsError &x)
1290 : {
1291 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
1292 : }
1293 0 : cfc.pa_p=Quantity(actualPA, "rad");
1294 : }
1295 0 : };
1296 :
1297 :
1298 : //
1299 : // Parser for parsing the NAME field of the SPECTRAL_WINDOW sub-table.
1300 : // Returns a vector of string containing tokens separated by "#"
1301 : // in the supplied band name string.
1302 : //
1303 0 : casacore::Vector<casacore::String> SynthesisUtils::parseBandName(const casacore::String& bandName)
1304 : {
1305 0 : casacore::Vector<casacore::String> tokens;
1306 :
1307 : // Allow a blank band name for older MSes where this is sometimes left blank.
1308 0 : if (bandName == "")
1309 : {
1310 0 : tokens.resize(1);tokens="";
1311 0 : return tokens;
1312 : }
1313 :
1314 : char *tok, *name;
1315 0 : int i=0;
1316 :
1317 0 : name = (char *)bandName.c_str();
1318 0 : if ((tok = std::strtok(name,"#"))!=NULL)
1319 : {
1320 0 : tokens.resize(i+1,true); tokens(i)="";
1321 0 : tokens(i++).assign(tok);
1322 : }
1323 :
1324 0 : while ((tok = std::strtok(NULL,"#"))!=NULL)
1325 : {
1326 0 : tokens.resize(i+1,true); tokens(i)="";
1327 0 : tokens(i++).assign(tok);
1328 : }
1329 0 : if (tokens(0)=="")
1330 : {
1331 0 : LogIO log_l(LogOrigin("SynthesisUtils", "parseBandName"));
1332 0 : log_l << "Error while parsing band name \"" << bandName << "\"" << LogIO::EXCEPTION;
1333 : }
1334 :
1335 : // for (i=0;i<3;i++)
1336 : // if (tokens(i)=="") log_l << "Error while parsing band name \"" << bandName << LogIO::EXCEPTION;
1337 0 : return tokens;
1338 : }
1339 :
1340 :
1341 : //
1342 : //-----------------------------------------------------------------------------------------
1343 : //
1344 0 : casacore::CoordinateSystem SynthesisUtils::makeModelGridFromImage(const std::string& modelImageName,
1345 : casacore::TempImage<casacore::DComplex>& modelImageGrid)
1346 : {
1347 : // This code is basically loading a floating point image from
1348 : // the disk, and loading it into a complex<double> image. This
1349 : // should be possible on-the-fly.
1350 : //
1351 : // However currently it is not possible to this OTF. So one has
1352 : // to load the image from the disk in a float image (copy-1 of
1353 : // the image in the memory). Then, since
1354 : // StokesImageUtil::From() does not work with complex<double>
1355 : // images, convert it first to a complex<float> image (equal to
1356 : // 2 more float buffers in the memory). And then covert the
1357 : // complex<float> image to a complex<Double> image (equal to 4
1358 : // more float buffers in the memory).
1359 : //
1360 : // In the end, just because of limitations of OTF type
1361 : // conversions, we end up with 7x memory footprint!
1362 :
1363 0 : casacore::LatticeBase* lattPtr = casacore::ImageOpener::openImage (modelImageName);
1364 : casacore::ImageInterface<float> *fImage;
1365 0 : fImage = dynamic_cast<casacore::ImageInterface<float>*>(lattPtr);
1366 :
1367 0 : TempImage<casacore::Complex> tmp(fImage->shape(), fImage->coordinates());
1368 0 : StokesImageUtil::From(tmp, *fImage);
1369 :
1370 0 : modelImageGrid = casacore::TempImage<casacore::DComplex> (fImage->shape(), fImage->coordinates());
1371 :
1372 : Bool d0,d1;
1373 0 : casacore::Array<casacore::DComplex> dcArray=modelImageGrid.get();
1374 0 : casacore::Array<casacore::Complex> fcArray=tmp.get();
1375 :
1376 0 : casacore::DComplex* dcArrayPtr= dcArray.getStorage(d0);
1377 0 : casacore::Complex* fcArrayPtr = fcArray.getStorage(d1);
1378 0 : IPosition ndx(4,0,0,0,0),shape=fImage->shape();
1379 :
1380 0 : for (ndx(0)=0; ndx(0)<shape(0);ndx(0)++)
1381 0 : for (ndx(1)=0; ndx(1)<shape(1);ndx(1)++)
1382 0 : for (ndx(2)=0; ndx(2)<shape(2);ndx(2)++)
1383 0 : for (ndx(3)=0; ndx(3)<shape(3);ndx(3)++)
1384 0 : dcArray(ndx) = DComplex(fcArray(ndx).real(), fcArray(ndx).imag());
1385 :
1386 0 : modelImageGrid.put(dcArray);
1387 0 : return fImage->coordinates();
1388 : }
1389 : //
1390 : //-------------------------------------------------------------------------------------
1391 : //
1392 0 : void SynthesisUtils::makeAWLists(const casacore::Vector<double>& wVals,
1393 : const casacore::Vector<double>& fVals,
1394 : const bool& wbAWP, const uint& nw,
1395 : const double& imRefFreq, const double& spwRefFreq,
1396 : casacore::Vector<int>& wNdxList,
1397 : casacore::Vector<int>& spwNdxList,
1398 : const int vbSPW)
1399 : {
1400 : //
1401 : // The following can be generalized to pick a subset of CFs along
1402 : // W and SPW axis in the CFB. Perhaps also useful in the longer
1403 : // run, e.g. when using a super-set CFC not all of which may be
1404 : // used for a given imaging.
1405 : //
1406 : // W-pixels in the CFC should be >= w-planes user setting
1407 0 : assert(wVals.nelements() >= nw);
1408 :
1409 : // Make list of W-CF indexes
1410 0 : int nWCFs=(nw<=1)?nw:wVals.nelements();
1411 0 : wNdxList.resize(nWCFs);
1412 0 : for(int i=0;i<nWCFs;i++) wNdxList[i] = i;
1413 :
1414 : // Make list of SPW-CF indexes
1415 0 : int nSPWCFs=fVals.nelements();
1416 0 : if (wbAWP==true)
1417 : {
1418 : // If a valid SPW ID is given, translate it to the spwNdx for the nearest SPW
1419 0 : if ((vbSPW>=0))// && (vbSPW <nSPWCFs))
1420 : {
1421 : int refSPW;
1422 0 : std::vector<double> stdList(nSPWCFs);
1423 0 : for (int i=0; i<nSPWCFs; i++) stdList[i] = fVals[i];
1424 : //Double refCFFreq =
1425 0 : SynthesisUtils::stdNearestValue(stdList, spwRefFreq, refSPW);
1426 :
1427 0 : spwNdxList.resize(1);
1428 0 : spwNdxList[0]=refSPW;
1429 : }
1430 : else
1431 : {
1432 0 : spwNdxList.resize(nSPWCFs);
1433 0 : for(int i=0;i<nSPWCFs;i++) spwNdxList[i] = i;
1434 : }
1435 : }
1436 : else
1437 : {
1438 : // For wbAWP=F, pick up the CF closest to the image reference frequency
1439 : int refSPW;
1440 0 : std::vector<double> stdList(nSPWCFs);
1441 0 : for (int i=0; i<nSPWCFs; i++) stdList[i] = fVals[i];
1442 : //Double refCFFreq =
1443 0 : SynthesisUtils::stdNearestValue(stdList, imRefFreq, refSPW);
1444 :
1445 0 : spwNdxList.resize(1);
1446 0 : spwNdxList[0]=refSPW;
1447 : }
1448 :
1449 0 : return;
1450 : }
1451 :
1452 : template
1453 : std::vector<Double>::iterator SynthesisUtils::Unique(std::vector<Double>::iterator first, std::vector<Double>::iterator last);
1454 : template
1455 : std::vector<Float>::iterator SynthesisUtils::Unique(std::vector<Float>::iterator first, std::vector<Float>::iterator last);
1456 : template
1457 : std::vector<Int>::iterator SynthesisUtils::Unique(std::vector<Int>::iterator first, std::vector<Int>::iterator last);
1458 :
1459 : template
1460 : Double SynthesisUtils::stdNearestValue(const vector<Double>& list, const Double& val, Int& index);
1461 : template
1462 : Float SynthesisUtils::stdNearestValue(const vector<Float>& list, const Float& val, Int& index);
1463 : template
1464 : Int SynthesisUtils::stdNearestValue(const vector<Int>& list, const Int& val, Int& index);
1465 :
1466 :
1467 : }
1468 :
1469 : //using namespace casacore;
1470 : } // namespace casa
|
