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