Line data Source code
1 : //# DelayRateFFT.cc: Part of the implementation of FringeJones
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 : #include <synthesis/MeasurementComponents/FringeJones.h>
28 :
29 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
30 : #include <casacore/lattices/Lattices/ArrayLattice.h>
31 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
32 :
33 : #include <casacore/casa/Exceptions/Error.h>
34 :
35 :
36 : #include <gsl/gsl_rng.h>
37 : #include <gsl/gsl_vector.h>
38 : #include <gsl/gsl_multimin.h>
39 :
40 : #include <iomanip> // needed for setprecision
41 :
42 : // DEVDEBUG gates the development debugging information to standard
43 : // error; it should be set to 0 for production.
44 :
45 : #define DEVDEBUG false
46 :
47 : using namespace casa::vi;
48 : using namespace casacore;
49 :
50 : namespace casa {
51 :
52 : Complex
53 : dotMatrixWithModel(const Matrix<Complex>& data, Double k, Double l, Double offset);
54 :
55 : Complex
56 : dotMatrixWithModel2(const Matrix<Complex>& f, Double k, Double l, Double offset);
57 :
58 : Complex
59 : dotWithOffsets(const Cube<Complex>& ft, const Vector<Float>& offsets, double k, double l);
60 :
61 : tuple<Double, Double, Double, Double>
62 : multibandFFTs(const Array<Complex>& ffts, const Vector<Float>& offsets);
63 :
64 : Double
65 : multibandFFT(const Vector<Complex>& peaks, const Vector<Float>& offsets);
66 :
67 : tuple<Double, Double, Double, Double>
68 : refineSearch2(const Cube<Complex>& ft, const Vector<Float>& offsets, Double ipkt0, Double ipkch0);
69 :
70 : tuple<Double, Double, Double, Double>
71 : bruteForceDelay(const Cube<Complex>& ft, const Vector<Float>& offsets, Int ipkt, Int ipkch);
72 :
73 :
74 :
75 125547 : static void unitize(Array<Complex>& vC)
76 : {
77 125547 : Array<Float> vCa(amplitude(vC));
78 : // Divide by non-zero amps
79 125547 : vCa(vCa<FLT_EPSILON)=1.0;
80 125547 : vC /= vCa;
81 125547 : }
82 :
83 8 : SDBListGridManagerCombo::SDBListGridManagerCombo(SDBList& sdbs) :
84 : SDBListGridManager(sdbs),
85 8 : nchan_(0)
86 : {
87 8 : std::set<Double> fmaxes_;
88 8 : std::set<Double> fmins_;
89 :
90 8 : if (sdbs.nSDB()==0) {
91 : // The for loop is fine with an empty list, but code below it
92 : // isn't and there's nothing to lose by bailing early!
93 : if (DEVDEBUG) {
94 : cerr << "No data I guess?";
95 : }
96 0 : return;
97 : }
98 :
99 800 : for (Int i=0; i != sdbs_.nSDB(); i++) {
100 792 : SolveDataBuffer& sdb = sdbs_(i);
101 792 : Int pspw = sdb.spectralWindow()(0);
102 792 : Double t = sdbs_(i).time()(0);
103 792 : times_.insert(t);
104 792 : if (spwins_.find(pspw) == spwins_.end()) {
105 8 : spwins_.insert(pspw);
106 8 : const Vector<Double>& fs = sdb.freqs();
107 8 : df_ = fs[1] - fs[0];
108 8 : pspwIdToFreqMap_[pspw] = &(sdb.freqs());
109 8 : nchan_ = max(nchan_, sdb.nChannels());
110 8 : fmaxes_.insert(fs(nchan_-1));
111 8 : fmins_.insert(fs(0));
112 : } else {
113 784 : continue;
114 : }
115 : }
116 8 : nt_ = times_.size();
117 : // nt_ = sdbs_.nSDB()/spwins_.size();
118 8 : tmin_ = *(times_.begin());
119 8 : tmax_ = *(times_.rbegin());
120 8 : dt_ = (tmax_ - tmin_)/(nt_ - 1);
121 8 : if (nchan_ == 1) {
122 0 : df_ = 1;
123 0 : return;
124 : }
125 8 : size_t i = 0;
126 16 : for (auto p=spwins_.begin(); p!=spwins_.end(); p++) {
127 8 : spwPMap_[*p] = i;
128 8 : i++;
129 : }
130 8 : }
131 :
132 :
133 :
134 :
135 : Float
136 208 : SDBListGridManagerCombo::getRefFreqFromLSPW(Int lspw) {
137 208 : auto p = spwins_.begin();
138 208 : std::advance(p, lspw);
139 208 : Int pspw = *p;
140 208 : const Vector<Double>& fs = *(pspwIdToFreqMap_[pspw]);
141 208 : return fs[0];
142 : }
143 :
144 225 : DelayRateFFT::DelayRateFFT(SDBList& sdbs, Int refant, Array<Double>& delayWindow, Array<Double>& rateWindow) :
145 450 : f0_(sdbs.centroidFreq() / 1.e9), // GHz, for delayrate calc
146 225 : refant_(refant),
147 225 : delayWindow_(delayWindow),
148 225 : rateWindow_(rateWindow),
149 225 : Vall_(),
150 225 : xcount_(),
151 225 : sumw_(),
152 225 : sumww_(),
153 225 : activeAntennas_(),
154 675 : allActiveAntennas_()
155 225 : {}
156 :
157 : DelayRateFFT *
158 225 : DelayRateFFT::makeAChild(int concat, SDBList& sdbs,
159 : Int refant,
160 : Array<Double>& delayWindow_,
161 : Array<Double>& rateWindow_)
162 : {
163 225 : if (concat) {
164 217 : return new DelayRateFFTConcat(sdbs, refant, delayWindow_, rateWindow_);
165 : } else {
166 8 : return new DelayRateFFTCombo(sdbs, refant, delayWindow_, rateWindow_);
167 : }
168 : }
169 :
170 : Matrix<Float>
171 0 : DelayRateFFT::delay() const {
172 0 : IPosition start(2, 1, 0);
173 0 : IPosition stop(2, 3*nCorr_-1, nElem_-1);
174 0 : IPosition stride(2, 3, 1);
175 0 : Slicer sl(start, stop, stride, Slicer::endIsLast);
176 0 : return param_(sl);
177 0 : }
178 :
179 : Matrix<Float>
180 0 : DelayRateFFT::rate() const {
181 0 : IPosition start(2, 2, 0);
182 0 : IPosition stop(2, 3*nCorr_-1, nElem_-1);
183 0 : IPosition stride(2, 3, 1);
184 0 : Slicer sl(start, stop, stride, Slicer::endIsLast);
185 0 : return param_(sl);
186 0 : }
187 :
188 :
189 362 : void DelayRateFFT::removeAntennasCorrelation(Int icor, std::set< Int > s) {
190 362 : std::set< Int > & as = activeAntennas_.find(icor)->second;
191 367 : for (std::set< Int >::iterator it=s.begin(); it!=s.end(); it++) {
192 5 : as.erase(*it);
193 : }
194 362 : }
195 :
196 : void
197 0 : DelayRateFFT::printActive() {
198 0 : for (Int icorr=0; icorr != nCorr_; icorr++) {
199 0 : cerr << "Antennas found for correlation " << icorr << ": ";
200 0 : for (std::set<Int>::iterator it = activeAntennas_[icorr].begin();
201 0 : it != activeAntennas_[icorr].end(); it++) {
202 0 : cerr << *it << ", ";
203 : }
204 0 : cerr << endl;
205 : }
206 0 : cerr << endl;
207 0 : }
208 :
209 :
210 :
211 : // DelayRateFFTCombo is modeled on DelayFFT in KJones.{cc|h}
212 8 : DelayRateFFTCombo::DelayRateFFTCombo(SDBList& sdbs, Int refant, Array<Double>& delayWindow,
213 8 : Array<Double>& rateWindow) :
214 : DelayRateFFT(sdbs, refant, delayWindow, rateWindow),
215 8 : gm_(sdbs)
216 : {
217 : if (DEVDEBUG) {
218 : gm_.describe();
219 : }
220 8 : nt_ = gm_.nt_;
221 8 : nChan_ = gm_.nchan_;
222 :
223 : // We might want to pad these too
224 8 : nPadFactor_ = 4;
225 8 : nPadT_ = nPadFactor_*nt_;
226 8 : nPadChan_ = nPadFactor_*nChan_;
227 8 : nspw_ = gm_.nSPW();
228 8 : dt_ = gm_.dt_;
229 8 : df_ = gm_.df_ / 1.e9;
230 :
231 8 : if (nt_ < 2) {
232 0 : throw(AipsError("Can't do a 2-dimensional FFT on a single timestep! Please consider changing solint to avoid orphan timesteps."));
233 : }
234 8 : IPosition ds = delayWindow_.shape();
235 8 : if (ds.size()!=1 || ds.nelements() != 1) {
236 0 : throw AipsError("delaywindow must be a list of length 2.");
237 : }
238 8 : IPosition rs = rateWindow_.shape();
239 8 : if (rs.size()!=1 || rs.nelements() != 1) {
240 0 : throw AipsError("ratewindow must be a list of length 2.");
241 : }
242 8 : Int nCorrOrig(sdbs(0).nCorrelations());
243 8 : nCorr_ = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
244 20 : for (Int i=0; i<nCorr_; i++) {
245 12 : activeAntennas_[i].insert(refant_);
246 : }
247 : // when we get the visCubecorrected it is already
248 : // reduced to parallel hands, but there isn't a
249 : // corresponding method for flags.
250 8 : Int corrStep = (nCorrOrig > 2 ? 3 : 1); // step for p-hands
251 800 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
252 792 : SolveDataBuffer& s(sdbs(ibuf));
253 36432 : for (Int irow=0; irow!=s.nRows(); irow++) {
254 35640 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
255 35640 : allActiveAntennas_.insert(a1);
256 35640 : allActiveAntennas_.insert(a2);
257 : }
258 : }
259 8 : nElem_ = 1 + *(allActiveAntennas_.rbegin()) ;
260 :
261 8 : IPosition aggregateDim(2, nCorr_, nElem_);
262 8 : xcount_.resize(aggregateDim);
263 8 : sumw_.resize(aggregateDim);
264 8 : sumww_.resize(aggregateDim);
265 8 : peak_.resize(aggregateDim);
266 :
267 8 : xcount_ = 0;
268 8 : sumw_ = 0.0;
269 8 : sumww_ = 0.0;
270 8 : IPosition dataSize(5, nCorr_, nElem_, nspw_, nPadT_, nPadChan_);
271 8 : Vall_.resize(dataSize);
272 8 : Int totalRows = 0;
273 8 : Int goodRows = 0;
274 :
275 : if (DEVDEBUG) {
276 : cerr << "DelayRateFFTCombo constructor: Here" << endl;
277 : }
278 800 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
279 792 : SolveDataBuffer& s(sdbs(ibuf));
280 792 : totalRows += s.nRows();
281 792 : if (!s.Ok())
282 0 : continue;
283 :
284 792 : Int nr = 0;
285 36432 : for (Int irow=0; irow!=s.nRows(); irow++) {
286 35640 : if (s.flagRow()(irow))
287 28512 : continue;
288 : Int iant;
289 35640 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
290 35640 : if (a1 == a2) { continue; } // we don't do autocorrelations
291 35640 : else if (a1 == refant_) { iant = a2; }
292 28512 : else if (a2 == refant_) { iant = a1; }
293 28512 : else { continue; } // not a baseline to reference antenna
294 : // v has shape (nelems, ?, nrows, nchannels); can't be a reference.
295 7128 : Cube<Complex> v = s.visCubeCorrected();
296 7128 : const Cube<Float>& w( s.weightSpectrum() );
297 7128 : const Cube<Bool>& fl( s.flagCube() );
298 7128 : Int pspw = s.spectralWindow()(0);
299 7128 : Int ispw = gm_.getLSPW(pspw);
300 7128 : Int t_index = gm_.getTimeIndex(s.time()(0));
301 : // FIXME! ispw is the Measurement Set spectral window; it
302 : // doesn't follow that that is addressible in our array
303 : // which is dimensioned for logical spectral
304 : // window. E.g. and i.e., if we are not joining spectral
305 : // windows, the spectral-window dimension of our target
306 : // array will always have size 1
307 7128 : IPosition start (5, 0, iant, ispw, t_index, 0);
308 7128 : IPosition stop (5, nCorr_, 1, 1, 1, nChan_);
309 7128 : IPosition stride(5, 1, 1, 1, 1, 1);
310 7128 : Slicer target_slice(start, stop, stride, Slicer::endIsLength);
311 : // Slicer::endIsLast is also possible
312 : if (DEVDEBUG && false) {
313 : cerr << " nSDBs" << sdbs.nSDB()
314 : << " nspwins " << gm_.spwins_.size() << endl;
315 : cerr << "nCorr_ " << nCorr_
316 : << " iant " << iant
317 : << " ispw " << ispw
318 : << " t_index " << t_index
319 : << " s.time()(0) " << s.time()(0)
320 : << " nt_ " << nt_
321 : << " dt_ " << dt_
322 : << " nChan_ " << nChan_
323 : << " corrStep " << corrStep
324 : << " distinct times " << gm_.times_.size()
325 : << endl;
326 : }
327 14256 : Slicer source_slice(IPosition(3, 0, 0, irow),
328 14256 : IPosition(3, nCorr_, nChan_, 1),
329 14256 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
330 :
331 14256 : Slicer flagSlice(IPosition(3, 0, 0, irow),
332 14256 : IPosition(3, nCorr_, nChan_, 1),
333 14256 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
334 7128 : nr++;
335 7128 : Array<Complex>rhs( v(source_slice).nonDegenerate(1) );
336 7128 : const Array<Float>& weights( w(source_slice).nonDegenerate(1) );
337 7128 : unitize(rhs);
338 7128 : Vall_(target_slice).nonDegenerate(1) = rhs * weights;
339 7128 : const Array<Bool>& flagged( fl(flagSlice).nonDegenerate(1) );
340 :
341 : // Zero flagged entries.
342 7128 : Vall_(target_slice).nonDegenerate(1)(flagged) = Complex(0.0);
343 :
344 7128 : if (!allTrue(flagged)) {
345 17226 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
346 : // Book keeping now done per spw!
347 10296 : IPosition p(3, icorr, iant, ispw);
348 10296 : Bool actually = false;
349 10296 : activeAntennas_[icorr].insert(iant);
350 586872 : for (Int ichan=0; ichan != (Int) nChan_; ichan++) {
351 576576 : IPosition pchan(2, icorr, ichan);
352 576576 : if (!flagged(pchan)) {
353 565488 : Float wv = weights(pchan);
354 565488 : xcount_(p)++;
355 565488 : sumw_(p) += wv;
356 565488 : sumww_(p) += wv*wv;
357 565488 : actually = true;
358 : }
359 576576 : }
360 10296 : if (actually) {
361 10098 : activeAntennas_[icorr].insert(iant);
362 10098 : goodRows++;
363 : }
364 10296 : }
365 : }
366 7128 : }
367 : }
368 8 : }
369 :
370 :
371 : // What even *is* this?
372 0 : DelayRateFFTCombo::DelayRateFFTCombo(Array<Complex>& data, Float f0, Float df, Float dt, SDBList& s,
373 0 : Array<Double>& delayWindow, Array<Double>& rateWindow) :
374 : DelayRateFFT(s, -1, delayWindow, rateWindow),
375 0 : gm_(s)
376 : {
377 0 : dt_ = dt;
378 0 : df_ = df;
379 0 : IPosition shape = data.shape();
380 0 : nCorr_ = shape(0);
381 0 : nElem_ = shape(1);
382 0 : nspw_ = shape(2);
383 0 : nt_ = shape(3);
384 0 : nChan_ = shape(4);
385 0 : IPosition dataSize(5, nCorr_, nElem_, nspw_, nt_, nChan_);
386 0 : Vall_.resize(dataSize);
387 :
388 0 : IPosition start(5, 0, 0, 0, 0, 0);
389 0 : IPosition stop(5, nCorr_, nElem_, nspw_, nt_, nChan_);
390 0 : IPosition stride(5, 1, 1, 1, 1);
391 0 : Slicer target_slice(start, stop, stride, Slicer::endIsLength);
392 0 : Vall_(target_slice) = data;
393 :
394 0 : unitize(Vall_);
395 :
396 0 : }
397 :
398 : void
399 8 : DelayRateFFTCombo::FFT() {
400 : if (DEVDEBUG) {
401 : cerr << "DelayRateFFTCombo::FFT()" << endl;
402 : }
403 : // Axes are 0: correlation (i.e., hand of polarization), 1: antenna, 2: spectral window, 3: time, 4: channel
404 : // the machinery is there to say that we only want to FFT the last two axes
405 8 : Vector<Bool> ax(5, false);
406 8 : ax(3) = true;
407 8 : ax(4) = true;
408 : // Also copied from DelayFFT in KJones.
409 : // we make a copy to FFT in place.
410 8 : ArrayLattice<Complex> c(Vall_);
411 8 : LatticeFFT::cfft0(c, ax, true);
412 : if (DEVDEBUG) {
413 : cerr << "FFT transformed" << endl;
414 : }
415 8 : }
416 :
417 : // In the new paradigm, this is where the new stuff happens. Instead of
418 : // interpolating the peaks on a single big grid, we have to synthesise
419 : // our estimate from separately FFT-ed spectral windows, using the new
420 : // off-grid peak formalism
421 : void
422 8 : DelayRateFFTCombo::searchPeak() {
423 : if (DEVDEBUG) {
424 : cerr << "DelayRateFFTCombo::searchPeak()" << endl;
425 : }
426 :
427 : // Recall param_ -> [phase, delay, rate] for each correlation
428 8 : param_.resize(3*nCorr_, nElem_); // This might be better done elsewhere.
429 8 : param_.set(0.0);
430 8 : flag_.resize(3*nCorr_, nElem_);
431 8 : flag_.set(true); // all flagged initially
432 :
433 8 : Double bw = Float(nChan_)*df_;
434 20 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
435 12 : flag_(icorr*3 + 0, refant()) = false;
436 12 : flag_(icorr*3 + 1, refant()) = false;
437 12 : flag_(icorr*3 + 2, refant()) = false;
438 132 : for (Int ielem=0; ielem<nElem_; ++ielem) {
439 120 : if (ielem==refant()) {
440 16 : continue;
441 : }
442 108 : if (activeAntennas_[icorr].find(ielem)==activeAntennas_[icorr].end()) {
443 4 : continue;
444 : }
445 :
446 : // Below is the gory details for turning delay window into index range
447 104 : Int sgn = (ielem < refant()) ? 1 : -1;
448 104 : Double d0 = sgn*delayWindow_(IPosition(1, 0));
449 104 : Double d1 = sgn*delayWindow_(IPosition(1, 1));
450 104 : if (d0 > d1) std::swap(d0, d1);
451 104 : d0 = max(d0, -0.5/df_);
452 104 : d1 = min(d1, (0.5-1/nChan_)/df_);
453 :
454 : // It's simpler to keep the ranges as signed integers and
455 : // handle the wrapping of the FFT in the loop over
456 : // indices. Recall that the FFT result returned has indices
457 : // that run from 0 to nChan_/2 -1 and then from
458 : // -nChan/2 to -1, so far as our delay is concerned.
459 104 : Int i0 = bw*d0;
460 104 : Int i1 = bw*d1;
461 104 : i0 *= nPadFactor_;
462 104 : i1 *= nPadFactor_;
463 104 : if (i1==i0) i1++;
464 : // Now for the gory details for turning rate window into index range
465 104 : Double width = nt_*dt_*1e9*f0_;
466 104 : Double r0 = sgn*rateWindow_(IPosition(1,0));
467 104 : Double r1 = sgn*rateWindow_(IPosition(1,1));
468 104 : if (r0 > r1) std::swap(r0, r1);
469 104 : r0 = max(r0, -0.5/(dt_*1e9*f0_));
470 104 : r1 = min(r1, (0.5 - 1/nt_)/(dt_*1e9*f0_));
471 :
472 104 : Int j0 = width*r0;
473 104 : Int j1 = width*r1;
474 104 : j0 *= nPadFactor_;
475 104 : j1 *= nPadFactor_;
476 104 : if (j1==j0) j1++;
477 : Array<Complex> blVis(
478 208 : Vall_(Slicer(
479 208 : IPosition(5, icorr, ielem, 0, 0, 0),
480 208 : IPosition(5, 1, 1, nspw_, nPadT_, nPadChan_),
481 104 : IPosition(5, 1, 1, 1, 1, 1),
482 208 : Slicer::endIsLength)).nonDegenerate(IPosition(1,2)));
483 :
484 104 : Vector<Float> offsets(nspw_);
485 104 : Double bw = gm_.nchan_ * gm_.df_;
486 208 : for (size_t lspw=0; lspw!=size_t(nspw_); lspw++) {
487 104 : Double f = gm_.getRefFreqFromLSPW(lspw);
488 104 : Double f0 = gm_.getRefFreqFromLSPW(0);
489 104 : offsets(lspw) = (f - f0)/bw;
490 : }
491 :
492 : // Get the DFT estimates of peak location instead
493 104 : tuple<Double, Double, Double, Double> mp = multibandFFTs(blVis , offsets);
494 104 : Double mpkt = std::get<0>(mp);
495 104 : Double mpkch = std::get<1>(mp);
496 104 : Double mpeak = std::get<2>(mp);
497 104 : Double marg = std::get<3>(mp);
498 :
499 104 : Complex c0 = dotWithOffsets(blVis, offsets, mpkt, mpkch);
500 : if (DEVDEBUG) {
501 : cerr << "DelayRateFFTCombo: abs(c1) " << abs(c0) << endl;
502 : }
503 104 : tuple<Double, Double, Double, Double> p = refineSearch(blVis, offsets, mpkt, mpkch);
504 104 : Double pkt = std::get<0>(p);
505 104 : Double pkch = std::get<1>(p);
506 104 : Double peak = std::get<2>(p);
507 104 : Double phase = std::get<3>(p);
508 :
509 : if (DEVDEBUG) {
510 : cerr << "[DelayRateFFTCombo::SearchPeak] " << "icorr " << icorr << " ielem " << ielem
511 : << " from (" << mpkt << ", " << mpkch << ", peak " << abs(c0)
512 : << ", angle " << arg(c0) << ")"
513 : << " to (" << pkt << ", " << pkch << ", peak " << peak << ", angle " << phase << ")"
514 : << endl;
515 : }
516 104 : peak_(IPosition(2, icorr, ielem)) = peak;
517 104 : param_(icorr*3 + 0, ielem) = sgn*phase;
518 104 : Double delay = (pkch)/Double(nChan_);
519 104 : if (delay > 0.5) delay -= 1.0; // fold
520 104 : delay /= df_; // nsec
521 104 : param_(icorr*3 + 1, ielem) = sgn*delay;
522 104 : Double rate = (pkt)/Double(nt_);
523 104 : if (rate > 0.5) rate -= 1.0;
524 104 : Double rate0 = rate/dt_;
525 104 : Double rate1 = rate0/(1e9 * f0_);
526 104 : param_(icorr*3 + 2, ielem) = Double(sgn*rate1);
527 : if (DEVDEBUG) {
528 : cerr << "delay " << sgn*delay << " rate1 " << sgn*rate1 << endl;
529 : }
530 : // Set 3 flags.
531 104 : flag_(icorr*3 + 0, ielem)=false;
532 104 : flag_(icorr*3 + 1, ielem)=false;
533 104 : flag_(icorr*3 + 2, ielem)=false;
534 : if (DEVDEBUG) {
535 : cerr << "End of this element/correlation combination" << endl;
536 : //cerr << "Set everything " << endl;
537 : }
538 104 : }
539 : }
540 8 : }
541 :
542 :
543 : tuple<Double, Double, Double, Double>
544 104 : multibandFFTs(const Array<Complex>& ffts, const Vector<Float>& offsets) {
545 : if (DEVDEBUG) {
546 : cerr << "multibandFFTs() " << endl;
547 : }
548 : // We take the individual band phases from the peaks of the SPWs
549 : // and assume they are the peculiar phases for their SPW; then we
550 : // calculate the Direct Discrete FT of those phases on their
551 : // offsets and read the peak off of that
552 :
553 : // This code is based on multibandFFT below, which we were using to
554 : // combine the peak values of individual SPW ffts; this version
555 : // combines *all* the values, because it turns out we can't
556 : // reliably get the peaks first.
557 :
558 104 : IPosition shp = ffts.shape();
559 104 : Int nspw = shp(0);
560 104 : Int nt = shp(1);
561 104 : Int nchan = shp(2);
562 :
563 104 : Int binScale = 1;
564 104 : size_t nbins= size_t(binScale*max(offsets)+1+0.5);
565 104 : Vector<Double> freqs(nbins);
566 104 : freqs = 0;
567 208 : for (size_t k=0; k!=nbins; k++) {
568 104 : freqs[k] = (Double(k) - nbins/2)/float(nbins);
569 : }
570 : // We need to store out results in an (nt, nchan, nbins) array this time
571 : // We do the loop manually because array arithmetic in Casacore is beyond me
572 104 : Array<Complex> X(IPosition(3, nbins, nt, nchan));
573 104 : X = 0;
574 :
575 208 : for (size_t ispw=0; ispw!=size_t(nspw); ispw++) {
576 41288 : for (size_t it=0; it != nt; it++) {
577 9266400 : for (size_t ichan=0; ichan != nchan; ichan++) {
578 9225216 : Double x = abs(ffts(IPosition(3, ispw, 0, 0)));
579 9225216 : if (std::isnan(x)) {cerr << "Skipping FFT " << ispw << endl;
580 0 : continue;
581 : }
582 18450432 : for (size_t ibin=0; ibin!=nbins; ibin++) {
583 9225216 : Complex rotation = exp(-Complex(0,1)*Complex((2.0*M_PI)*offsets(ispw)*freqs(ibin)));
584 9225216 : X(IPosition(3, ibin, it, ichan)) += ffts(IPosition(3, ispw, it, ichan))*rotation;
585 : }
586 : }
587 : }
588 : }
589 : /*
590 : // New!
591 : Array<Complex> rateAveraged(IPosition(2, nbins, nchan));
592 : rateAveraged = 0;
593 : for (size_t ibin; ibin!=nbins; ibin++) {
594 : for (size_t ichan=0; ichan != nchan; ichan++) {
595 : for (size_t it=0; it != nt; it++) {
596 : rateAveraged(IPosition(2, ibin, ichan)) += abs(ffts(IPosition(3, ibin, it, ichan)));
597 : }
598 : }
599 : }
600 :
601 : // Search for a peak with delay in the average first
602 : Int ichanmax1 = -1;
603 : Int ibinmax1 = -1;
604 : Double zmax1 = -1.0;
605 : for (Int ichan=0; ichan != nchan; ichan++) {
606 : for (Int ibin=0; ibin != nbins; ibin++) {
607 : Complex c = rateAveraged(IPosition(2, ibin, ichan));
608 : Double a = abs(c);
609 : if (a>zmax1) {
610 : ichanmax1 = ichan;
611 : ibinmax1 = ibin;
612 : zmax1 = a;
613 : }
614 : }
615 : }
616 : // Then we search the time dimension
617 : zmax1 = -1.0;
618 : Int itmax1 = -1;
619 : for (Int it = 0; it != nt; it++) {
620 : Complex c = X(IPosition(3, ibinmax1, it, ichanmax1));
621 : Double a = abs(c);
622 : if (a>zmax1) {
623 : itmax1 = it;
624 : }
625 : }
626 : cerr << "Averaged stacked DFT peaks at " << itmax1 << "/" << nt << ", "
627 : << ichanmax1 << "/" << nchan << "; Bin " << ibinmax1
628 : << " Peak " << abs(X(IPosition(3, ibinmax1, itmax1, ichanmax1)))
629 : << endl;
630 : */
631 : // We search for the maximum ourself, by brute force.
632 104 : Int itmax = -1;
633 104 : Int ichanmax = -1;
634 104 : Int ibinmax = -1;
635 104 : Double zmax = -1.0;
636 208 : for (Int ibin=0; ibin != nbins; ibin++) {
637 41288 : for (Int it = 0; it != nt; it++) {
638 9266400 : for (Int ichan=0; ichan != nchan; ichan++) {
639 9225216 : Complex c = X(IPosition(3, ibin, it, ichan));
640 9225216 : Double a = abs(c);
641 9225216 : if (a>zmax) {
642 116 : itmax = it;
643 116 : ichanmax = ichan;
644 116 : ibinmax = ibin;
645 116 : zmax = a;
646 : }
647 : }
648 : }
649 : }
650 : //cerr << "Unaveraged stacked DFT peaks at " << itmax << "/" << nt << ", "
651 : // << ichanmax << "/" << nchan << "; Bin " << ibinmax
652 : // << " Peak " << abs(X(IPosition(3, ibinmax, itmax, ichanmax)))
653 : // << endl;
654 : //for (Int i=0; i < nspw; i++) {
655 : // cerr << "Peak for bin " << i << " = " << abs(ffts(IPosition(3, i, itmax, ichanmax))) << endl;
656 : //}
657 :
658 104 : Complex c = X(IPosition(3, ibinmax, itmax, ichanmax));
659 104 : Double dpkch = freqs(ibinmax);
660 :
661 104 : tuple<Double, Double, Double, Double> t;
662 104 : t = std::make_tuple(Double(itmax), Double(ichanmax) + dpkch, abs(c), arg(c));
663 : // t = std::make_tuple(Double(itmax), Double(ichanmax), abs(c), arg(c));
664 : if (DEVDEBUG) {
665 : cerr << "itmax " << itmax << " ichanmax " << ichanmax << " ibinmax " << ibinmax
666 : << " peak " << abs(c) << endl;
667 : cerr << "MultibandFFTs dpkch=" << dpkch << endl;
668 : cerr << "multibandFFTs() finished" << endl;
669 : }
670 208 : return t;
671 104 : }
672 :
673 : Double
674 0 : multibandFFT(const Vector<Complex>& peaks, const Vector<Float>& offsets) {
675 : // We take the individual band phases from the peaks of the SPWs
676 : // and assume they are the peculiar phases for their SPW; then we
677 : // calculate the Direct Discrete FT of those phases on their
678 : // offsets and read the peak off of that
679 0 : Int nspw = peaks.size();
680 :
681 0 : Vector<Float> amps( amplitude(peaks) );
682 : Int len;
683 0 : amps.shape(len);
684 :
685 0 : Float eps = 1e-10;
686 0 : std::set<size_t> badInds;
687 0 : for (size_t i=0; i!=size_t(len); i++) {
688 0 : if (amps[i] < eps) {
689 : if (DEVDEBUG) {
690 : cerr << " blacklisting " << i << endl;
691 : }
692 0 : badInds.insert(i);
693 : }
694 : }
695 0 : Vector<Complex> phasors( peaks/amplitude(peaks) );
696 :
697 : // I want to print angles but I can't get the API to do that
698 : // in a reasonable way; arg(phasors) doesn't work
699 : // Vector<Float> angles ( arg(phasors) );
700 : if (DEVDEBUG) {
701 : cerr << "phasors " << phasors << endl;
702 : }
703 : // FIXME: the last offset is at the *beginning* of the subband!
704 : // size_t nbins=2*size_t(max(offsets)+1+0.5);
705 0 : size_t nbins= size_t(max(offsets)+1+0.5);
706 :
707 0 : Vector<Float> freqs(nbins);
708 0 : freqs = 0;
709 0 : for (size_t k=0; k!=nbins; k++) {
710 0 : freqs[k] = (Float(k) - nbins/2)/float(nbins);
711 : }
712 0 : Vector<Complex> X(nbins);
713 0 : X = 0;
714 0 : for (size_t n=0; n!=size_t(nspw); n++) {
715 0 : if (badInds.find(n) != badInds.end()) {
716 : if (DEVDEBUG) {
717 : cerr << " skipping " << n << endl;
718 : }
719 0 : continue;
720 : }
721 0 : for (size_t k=0; k!=nbins; k++) {
722 : // FIXME: also swap signs here!
723 0 : X[k] += peaks[n]*exp(-Complex(0,1)*Complex((2.0*M_PI)*offsets(n)*freqs(k)));
724 : // X[k] += phasors[n]*exp(-Complex(0,1)*Complex((2.0*M_PI)*offsets(n)*freqs(k)));
725 : }
726 : }
727 0 : Int k_max = -1;
728 0 : Float zmax = -1.0;
729 0 : for (size_t k=0; k!=nbins; k++) {
730 0 : Complex z = X[k];
731 0 : Float a = abs(z);
732 : if (DEVDEBUG) {
733 : cerr << "DFT Bin " << k << " freq " << freqs(k) << " coeff " << a << endl;
734 : }
735 0 : if (a>zmax) {
736 0 : k_max = k;
737 0 : zmax = a;
738 : }
739 : }
740 : if (DEVDEBUG) {
741 : cerr << "k_max = " << k_max
742 : << "; peak at k_max= " << abs(X[k_max])
743 : << "; freqs(k_max) = " << freqs(k_max) << endl;
744 : }
745 0 : return freqs(k_max);
746 0 : }
747 :
748 : // Auxilliary function
749 : Complex
750 1644 : dotWithOffsets(const Cube<Complex>& ft, const Vector<Float>& offsets, double k, double l) {
751 1644 : Int nspw = ft.nrow();
752 1644 : Int ni = ft.ncolumn();
753 1644 : Int nj = ft.nplane();
754 :
755 1644 : if (k<0) k += (ni);
756 1644 : if (l<0) l += (nj);
757 1644 : Complex p(0.0, 0.0);
758 3288 : for (Int s=0; s!=nspw; s++) {
759 1644 : const Matrix<Complex>& ft_s = ft.yzPlane(s);
760 1644 : Double f_off = offsets(s);
761 1644 : Complex r = dotMatrixWithModel(ft_s, k, l, f_off);
762 1644 : p += r;
763 1644 : }
764 1644 : return p;
765 : }
766 :
767 : // We need a function to minimize, which has to return a real value,
768 : // but when we're done we need to find the peak and also its argument,
769 : // so we factor out the complex version...
770 : Complex
771 1540 : c_peak_fn(const gsl_vector *x, void *vparams) {
772 1540 : std::pair< Cube<Complex> const *, Vector<Float> const * > *params =
773 : (std::pair< Cube<Complex> const *, Vector<Float> const * > *)(vparams );
774 :
775 1540 : double k = gsl_vector_get(x, 0);
776 1540 : double l = gsl_vector_get(x, 1);
777 :
778 1540 : const Cube<Complex>& ft ( *(params->first ));
779 1540 : const Vector<Float>& offsets ( *(params->second));
780 1540 : return dotWithOffsets(ft, offsets, k, l);
781 : }
782 :
783 : // ... and then call it from a real version that we can give to the
784 : // optimizer (with a sign change; by tradition these are minimizer
785 : // routines and we want a maximum)
786 : double
787 1436 : my_peak_fn(const gsl_vector *x, void *vparams) {
788 1436 : Complex p = c_peak_fn(x, vparams);
789 1436 : return -abs(p);
790 : }
791 :
792 : tuple<Double, Double, Double, Double>
793 0 : bruteForceDelay(const Cube<Complex>& ft, const Vector<Float>& offsets, Int ipkt, Int ipkch) {
794 0 : Double k(ipkt);
795 0 : Double l(ipkch);
796 0 : Int n = 10;
797 :
798 0 : Double l_p = l;
799 0 : Complex peak = Complex(0, 0);
800 0 : for (Int i=-n; i!=n+1; i++) {
801 0 : Double dch = double(i)/(2*n);
802 0 : Complex p = dotWithOffsets(ft, offsets, k, l+dch);
803 0 : if (abs(p) > abs(peak)) {
804 : if (DEVDEBUG) {
805 : cerr << "[bruteForceDelay] l " << l+dch << " peak " << abs(p) << endl;
806 : }
807 0 : peak = p;
808 0 : l_p = l+dch;
809 :
810 : }
811 : }
812 0 : tuple<Double, Double, Double, Double> t;
813 0 : t = std::make_tuple(k, l_p, abs(peak), arg(peak));
814 0 : return t;
815 : }
816 :
817 : tuple<Double, Double, Double, Double>
818 104 : DelayRateFFTCombo::refineSearch(const Cube<Complex>& ft, const Vector<Float>& offsets, Double pkt, Double pkch) {
819 : // small@jive.eu: I borrowed most of this code from the GSL documentation of multimin:
820 : // <https://www.gnu.org/software/gsl/doc/html/multimin.html>
821 : if (DEVDEBUG) {
822 : cerr << "refineSearch" << endl;
823 : }
824 104 : const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
825 : /* Starting point */
826 104 : gsl_vector *x = gsl_vector_alloc (2);
827 104 : gsl_vector_set(x, 0, pkt);
828 104 : gsl_vector_set(x, 1, pkch);
829 : if (DEVDEBUG) {
830 : cerr << "pkch starts at " << pkch << endl;
831 : }
832 : /* Set initial step sizes */
833 : /* Fixme! I need to think harder about this value! */
834 104 : gsl_vector* steps = gsl_vector_alloc (2);
835 104 : gsl_vector_set_all(steps, 0.01);
836 :
837 : /* Initialize method and iterate */
838 :
839 : gsl_multimin_function minex_func;
840 104 : std::pair< Cube<Complex> const * , Vector<Float> const * > par = make_pair(&ft, &offsets);
841 :
842 104 : minex_func.n = 2;
843 104 : minex_func.f = my_peak_fn;
844 104 : minex_func.params = ∥
845 :
846 104 : gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, 2);
847 104 : gsl_multimin_fminimizer_set (s, &minex_func, x, steps);
848 :
849 104 : float minstep = 1e-5;
850 : int status;
851 104 : size_t iter = 0;
852 : do {
853 : if (DEVDEBUG) {
854 : cerr << "Starting iteration" << endl;
855 : }
856 :
857 520 : iter++;
858 520 : status = gsl_multimin_fminimizer_iterate(s);
859 520 : if (status) break;
860 520 : double size = gsl_multimin_fminimizer_size(s);
861 520 : status = gsl_multimin_test_size(size, minstep);
862 520 : if (status == GSL_SUCCESS) {
863 0 : printf ("converged to minimum at\n");
864 : }
865 : if (DEVDEBUG) {
866 : printf("%5zd ipkt %12.5e ipkch %12.5e f() = %12.5g size = %12.5f\n",
867 : iter,
868 : gsl_vector_get(s->x, 0),
869 : gsl_vector_get(s->x, 1),
870 : s->fval,
871 : size);
872 : }
873 520 : } while (status == GSL_CONTINUE && iter < 5);
874 : // while (status == GSL_CONTINUE && iter < 100);
875 104 : tuple<Double, Double, Double, Double> p;
876 : // if (status == GSL_SUCCESS) {
877 : if (1) {
878 : // Double peak = gsl_multimin_fminimizer_minimum(s);
879 104 : Double ipkt = gsl_vector_get(s->x, 0);
880 104 : Double ipkch = gsl_vector_get(s->x, 1);
881 104 : Complex c = c_peak_fn(s->x, &par);
882 104 : p = std::make_tuple(ipkt, ipkch, abs(c), arg(c));
883 : } else {
884 : p = std::make_tuple(pkt, pkch, 0.0, 0.0);
885 : }
886 104 : gsl_vector_free(steps);
887 104 : gsl_vector_free(x);
888 104 : gsl_multimin_fminimizer_free(s);
889 :
890 208 : return p;
891 : }
892 :
893 : Float
894 104 : DelayRateFFTCombo::snr(Int icorr, Int ielem, Float delay, Float rate) {
895 : // We calculate a signal-to-noise ration for the 2D FFT fringefit
896 : // using a formula transcribed from AIPS FRING.
897 : //
898 : // Have to convert delay and rate back into indices on the padded 2D grid.
899 :
900 104 : IPosition p(2, icorr, ielem);
901 104 : Float peak = peak_(p);
902 104 : if (peak > 0.999*sumw_(p)) {
903 : //cerr << "Clipping peak for element " << ielem << " from " << peak << " to " << 0.999*sumw_(p) << endl;
904 24 : peak=0.999*sumw_(p);
905 : }
906 : // xcount is number of data points for baseline to ielem
907 : // sumw is sum of weights,
908 : // sumww is sum of squares of weights
909 : Float cwt;
910 104 : if (fabs(sumw_(p))<FLT_EPSILON) {
911 4 : cwt = 0;
912 : if (DEVDEBUG) {
913 : cerr << "Correlation " << icorr << " antenna " << ielem << ": sum of weights is zero." << endl;
914 : }
915 : } else {
916 100 : Float x = M_PI/2*peak/sumw_(p);
917 : // The magic numbers in the following formula are from AIPS FRING
918 100 : cwt = (pow(tan(x), 1.163) * sqrt(sumw_(p)/sqrt(sumww_(p)/xcount_(p))));
919 : if (DEVDEBUG) {
920 : cerr << "Correlation " << icorr << " antenna " << ielem
921 : << " peak=" << peak << "; xang=" << x << "; xcount=" << xcount_(p) << "; sumw=" << sumw_(p) << "; sumww=" << sumww_(p)
922 : << " snr " << cwt << endl;
923 : }
924 : }
925 104 : return cwt;
926 104 : }
927 :
928 :
929 : // There isn't a usable sinc lying around that I can see
930 0 : Double sinc(Double x) {
931 0 : Double p = M_PI*x;
932 0 : return sin(p)/p;
933 : }
934 :
935 : // This implements a Fancy Sinc strategy for finding the peak of a
936 : // single SPW. I no longer know if this can be extended to multiple
937 : // SPWs.
938 : tuple<Double, Double, Double, Double>
939 0 : refineSearch2(const Cube<Complex>& ft, const Vector<Float>& offsets, Int ipkt0, Int ipkch0) {
940 : Double imax, jmax;
941 : // FIXME! We need to do a single SPW first
942 0 : if (ft.nrow() > 1) {
943 0 : throw AipsError("Only one spw for now");
944 : }
945 0 : Int s = 0;
946 0 : const Matrix<Complex>& ft_s(ft.yzPlane(s));
947 0 : size_t ni = ft_s.nrow();
948 0 : size_t nj = ft_s.ncolumn();
949 : // we need to do a roll left/down and sum but Casacore doesn't have a matrix roll so we may as well do it by hand
950 0 : Float pmax = -1;
951 0 : size_t ipkt(0);
952 0 : size_t ipkch(0);
953 0 : Matrix<Float> sumAbs(ni, nj);
954 0 : for (size_t i=0; i!=ni; i++) {
955 : // FIXME: The variable i1 is unused!
956 : // Either use it or lose it!
957 : // size_t i1 = (i==ni-1) ? 0 : i+1;
958 0 : for (size_t j=0; j!=nj; j++) {
959 0 : size_t j1 = (j==nj-1) ? 0 : j+1;
960 0 : Float sumAbs = (abs(ft_s(i, j)) + abs(ft_s(i, j1)));
961 0 : if ((sumAbs) > pmax) {
962 0 : ipkt = i;
963 0 : ipkch = j;
964 0 : pmax = sumAbs;
965 : }
966 : }
967 : }
968 0 : Double y0 = abs(ft_s(ipkt, ipkch));
969 0 : Double y1 = abs(ft_s(ipkt, (ipkch==nj-1)? 0 : ipkch+1));
970 0 : Double ym1 = abs(ft_s(ipkt, (ipkch==0)? nj-1 : ipkch-1));
971 :
972 0 : Double z1 = abs(ft_s((ipkt==ni-1) ? 0 : ipkt+1, ipkch));
973 0 : Double zm1 = abs(ft_s((ipkt==0) ? ni-1 : ipkt-1, ipkch));
974 :
975 :
976 0 : Double d0 = y1/(y0+y1);
977 0 : Double peak( y0/sinc(d0) );
978 : if (DEVDEBUG) {
979 : cerr << "y0 " << y0 << " y1 " << y1 << " ym1 " << ym1 << " d0 " << d0 << " sinc(d0) "
980 : << sinc(d0) << " peak " << peak << endl;
981 : cerr << "y0 " << y0 << " z1 " << z1 << " zm1 " << zm1 << endl;
982 : cerr << real(abs(ft_s.row(ipkt)))/y0 << endl;
983 : }
984 0 : Double phase( 0.0 );
985 0 : imax = Double(ipkt);
986 0 : jmax = Double(ipkch) + d0;
987 0 : tuple<Double, Double, Double, Double> p = std::make_tuple(imax, jmax, peak, phase);
988 0 : return p;
989 0 : }
990 :
991 : Complex
992 1644 : dotMatrixWithModel(const Matrix<Complex>& data, Double k, Double l, Double offset)
993 : {
994 1644 : size_t ni = data.nrow();
995 1644 : size_t nj = data.ncolumn();
996 1644 : Matrix<Complex> model(ni, nj);
997 :
998 : // Note that k is the time index of the array, l is the frequency index.
999 : // Offsets correspond to frequencies
1000 1644 : Double eps = 1e-8;
1001 1644 : Int k_int = floor(k);
1002 1644 : Double k_del = k - k_int;
1003 1644 : Bool k_flag = (fabs(k_del) < eps);
1004 1644 : Int l_int = floor(l);
1005 1644 : Double l_del = l - l_int;
1006 1644 : Bool l_flag = (fabs(l_del) < eps);
1007 1644 : Complex t0;
1008 1644 : Complex t1;
1009 : // FIXME! We're messing with reversing signs
1010 652668 : for (size_t i=0; i!=ni; i++) {
1011 651024 : if (k_flag) {
1012 127512 : t0 = Complex(i==k);
1013 : } else { // if k isn't an integer!
1014 523512 : t0 = ( (1-exp(Complex(0, -1.0*(2.0*M_PI)*(k-i))))/
1015 1047024 : (1-exp(Complex(0, -1.0*(2.0*M_PI)*(k-i)/Double(ni)))) );
1016 523512 : t0 /= ni;
1017 : }
1018 146480400 : for (size_t j=0; j!=nj; j++) {
1019 145829376 : if (l_flag) {
1020 28562688 : t1 = Complex(j==l);
1021 : } else { // if l isn't an integer!
1022 117266688 : t1 = ( (1-exp(Complex(0, -1.0*(2.0*M_PI)*(l-j))))/
1023 234533376 : (1-exp(Complex(0, -1.0*(2.0*M_PI)*(l-j)/Double(nj)))) );
1024 117266688 : t1 /= nj;
1025 : }
1026 145829376 : model(i, j) = exp(Complex(0, -1.0*(2.0*M_PI)*offset*l_del))*t0*t1;
1027 : }
1028 : }
1029 1644 : Complex t2 = sum(data*model);
1030 1644 : return t2;
1031 1644 : }
1032 :
1033 :
1034 : Complex
1035 0 : dotMatrixWithModel2(const Matrix<Complex>& data, Double k0, Double l0, Double offset)
1036 : {
1037 0 : size_t ni = data.nrow();
1038 0 : size_t nj = data.ncolumn();
1039 0 : Matrix<Complex> model(ni, nj);
1040 0 : model = 0.0;
1041 :
1042 0 : Double eps = 1e-8;
1043 0 : Int k_int = Int(ceil(k0));
1044 0 : Int l_int = Int(ceil(l0));
1045 :
1046 0 : Double d0 = k_int - k0;
1047 0 : Double d1 = l_int - l0;
1048 :
1049 0 : Bool k_flag = (fabs(d0) < eps);
1050 0 : Bool l_flag = (fabs(d1) < eps);
1051 :
1052 0 : Complex c0, c1;
1053 :
1054 0 : Complex s (0.0);
1055 0 : Complex t0 = Complex(sin(M_PI*d0)/M_PI)*exp(Complex(0, M_PI*d0));
1056 0 : Complex t1 = Complex(sin(M_PI*d1)/M_PI)*exp(Complex(0, M_PI*d1));
1057 0 : for (Int dk=-2; dk!=2; dk++) {
1058 0 : size_t k = (k_int + dk + ni) % ni;
1059 0 : if (k_flag) {
1060 0 : c0 = Complex(dk==0);
1061 : } else { // if k isn't an integer!
1062 0 : c0 = Complex(1/(d0+dk))*t0;
1063 : }
1064 0 : for (Int dl=-2; dl!=2; dl++) {
1065 0 : size_t l = (l_int+dl+nj) % nj;
1066 0 : if (l_flag) {
1067 0 : c1 = Complex(dl==0);
1068 : } else { // if l isn't an integer!
1069 0 : c1 = Complex(1/(d1+dl))*t1;
1070 : }
1071 0 : Complex d = data(k, l);
1072 0 : Complex m = c0*c1;
1073 0 : s += d*m;
1074 : }
1075 : }
1076 0 : return s;
1077 0 : }
1078 :
1079 217 : SDBListGridManagerConcat::SDBListGridManagerConcat(SDBList& sdbs) :
1080 217 : SDBListGridManager(sdbs)
1081 : {
1082 217 : std::set<Double> fmaxes;
1083 217 : std::set<Double> fmins;
1084 217 : Float dfn(0.0);
1085 217 : Int totalChans0(0);
1086 :
1087 217 : if (sdbs_.nSDB()==0) {
1088 : // The for loop is fine with an empty list, but code below it
1089 : // isn't and there's nothing to lose by bailing early!
1090 0 : return;
1091 : }
1092 :
1093 26944 : for (Int i=0; i != sdbs_.nSDB(); i++) {
1094 26727 : SolveDataBuffer& sdb = sdbs_(i);
1095 :
1096 26727 : Int spw = sdb.spectralWindow()(0);
1097 26727 : Double t = sdbs_(i).time()(0);
1098 26727 : times_.insert(t);
1099 : // If we haven't seen this spectral window before, handle it
1100 26727 : if (spwins_.find(spw) == spwins_.end()) {
1101 229 : spwins_.insert(spw);
1102 229 : Int nActualChannels = sdb.nChannels();
1103 229 : setNChans(spw, nActualChannels);
1104 229 : const Vector<Double>& fs = sdb.freqs();
1105 229 : spwIdToFreqMap_[spw] = &(sdb.freqs());
1106 229 : fmaxes.insert(fs(nActualChannels-1));
1107 229 : fmins.insert(fs(0));
1108 229 : totalChans0 += nActualChannels;
1109 229 : dfn = (fs(nActualChannels-1) - fs(0))/(nActualChannels-1);
1110 : } else {
1111 26498 : continue;
1112 : }
1113 : }
1114 217 : nt_ = times_.size();
1115 217 : tmin_ = *(times_.begin());
1116 217 : tmax_ = *(times_.rbegin());
1117 217 : dt_ = (tmax_ - tmin_)/(nt_ - 1);
1118 217 : fmin_ = *(fmins.begin());
1119 217 : fmax_ = *(fmaxes.rbegin());
1120 217 : totalChans_ = round((fmax_ - fmin_)/dfn + 1);
1121 217 : df_ = (fmax_ - fmin_)/(totalChans_-1);
1122 : if (DEVDEBUG) {
1123 : cerr << "Global fmin " << fmin_ << " global max " << fmax_ << endl;
1124 : cerr << "nt " << nt_ << " dt " << dt_ << endl;
1125 : cerr << "tmin " << tmin_ << " tmax " << tmax_ << endl;
1126 : cerr << "Global dt " << tmax_ - tmin_ << endl;
1127 : cerr << "Global df " << df_ << endl;
1128 : cerr << "I guess we'll need " << totalChans_ << " freq points in total." << endl;
1129 : cerr << "Compared to " << totalChans0 << " with simple-minded concatenation." << endl;
1130 : cerr << "spwins_.size() " << spwins_.size() << endl;
1131 : cerr << "nSPW() " << nSPW() << endl;
1132 : }
1133 217 : }
1134 :
1135 : // checkAllGridPoints is a diagnostic funtion that should not be called
1136 : // in production releases, but it doesn't do any harm to have it
1137 : // latent.
1138 : void
1139 0 : SDBListGridManagerConcat::checkAllGridpoints() {
1140 0 : map<Int, Vector<Double> const *>::iterator it;
1141 0 : for (it = spwIdToFreqMap_.begin(); it != spwIdToFreqMap_.end(); it++) {
1142 0 : Int spwid = it->first;
1143 0 : Vector<Double> const* fs = it->second;
1144 : Int length;
1145 0 : fs->shape(length);
1146 0 : for (Int i=0; i!=length; i++) {
1147 0 : Double f = (*fs)(i);
1148 0 : Int j = bigFreqGridIndex(f);
1149 : if (DEVDEBUG) {
1150 : cerr << "spwid, i = (" << spwid << ", " << i << ") => " << j << " (" << f << ")" << endl;
1151 : }
1152 : }
1153 : }
1154 : if (DEVDEBUG) {
1155 : cerr << "[1] spwins.size() " << nSPW() << endl;
1156 : }
1157 0 : }
1158 :
1159 : Int
1160 118419 : SDBListGridManagerConcat::swStartIndex(Int spw) {
1161 118419 : Vector<Double> const* fs = spwIdToFreqMap_[spw];
1162 118419 : Double f0 = (*fs)(0);
1163 118419 : return bigFreqGridIndex(f0);
1164 : }
1165 :
1166 217 : DelayRateFFTConcat::DelayRateFFTConcat(SDBList& sdbs, Int refant, Array<Double>& delayWindow,
1167 217 : Array<Double>& rateWindow) :
1168 : DelayRateFFT(sdbs, refant, delayWindow, rateWindow),
1169 217 : gm_(sdbs)
1170 : {
1171 : if (DEVDEBUG) {
1172 : gm_.describe();
1173 : cerr << "gm_.nSPW() " << gm_.nSPW() << endl;
1174 : }
1175 217 : nPadFactor_= max(2, 8 / gm_.nSPW());
1176 217 : nt_ = gm_.nt_;
1177 217 : nPadT_ = nPadFactor_ * nt_;
1178 217 : nChan_ = gm_.nChannels();
1179 217 : nPadChan_ = nPadFactor_*nChan_;
1180 217 : dt_ = gm_.dt_;
1181 217 : f0_ = sdbs.centroidFreq() / 1.e9; // GHz, for delayrate calc
1182 217 : df_ = gm_.df_ / 1.e9;
1183 217 : df_all_ = gm_.fmax_ - gm_.fmin_;
1184 : // This check should be commented out in production:
1185 : // gm_.checkAllGridpoints();
1186 217 : if (nt_ < 2) {
1187 0 : throw(AipsError("Can't do a 2-dimensional FFT on a single timestep! Please consider changing solint to avoid orphan timesteps."));
1188 : }
1189 217 : IPosition ds = delayWindow_.shape();
1190 217 : if (ds.size()!=1 || ds.nelements() != 1) {
1191 0 : throw AipsError("delaywindow must be a list of length 2.");
1192 : }
1193 217 : IPosition rs = rateWindow_.shape();
1194 217 : if (rs.size()!=1 || rs.nelements() != 1) {
1195 0 : throw AipsError("ratewindow must be a list of length 2.");
1196 : }
1197 217 : Int nCorrOrig(sdbs(0).nCorrelations());
1198 217 : nCorr_ = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
1199 567 : for (Int i=0; i<nCorr_; i++) {
1200 350 : activeAntennas_[i].insert(refant_);
1201 : }
1202 :
1203 : // when we get the visCubecorrected it is already
1204 : // reduced to parallel hands, but there isn't a
1205 : // corresponding method for flags.
1206 217 : Int corrStep = (nCorrOrig > 2 ? 3 : 1); // step for p-hands
1207 26944 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
1208 26727 : SolveDataBuffer& s(sdbs(ibuf));
1209 527557 : for (Int irow=0; irow!=s.nRows(); irow++) {
1210 500830 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
1211 500830 : allActiveAntennas_.insert(a1);
1212 500830 : allActiveAntennas_.insert(a2);
1213 : }
1214 : }
1215 217 : nElem_ = 1 + *(allActiveAntennas_.rbegin()) ;
1216 217 : IPosition aggregateDim(2, nCorr_, nElem_);
1217 217 : xcount_.resize(aggregateDim);
1218 217 : sumw_.resize(aggregateDim);
1219 217 : sumww_.resize(aggregateDim);
1220 :
1221 217 : xcount_ = 0;
1222 217 : sumw_ = 0.0;
1223 217 : sumww_ = 0.0;
1224 :
1225 : if (DEVDEBUG) {
1226 : cerr << "Filling FFT grid with " << sdbs.nSDB() << " data buffers." << endl;
1227 : }
1228 217 : IPosition paddedDataSize(4, nCorr_, nElem_, nPadT_, nPadChan_);
1229 217 : Vpad_.resize(paddedDataSize);
1230 217 : Int totalRows = 0;
1231 217 : Int goodRows = 0;
1232 :
1233 :
1234 26944 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
1235 26727 : SolveDataBuffer& s(sdbs(ibuf));
1236 26727 : totalRows += s.nRows();
1237 26727 : if (!s.Ok())
1238 30 : continue;
1239 :
1240 26697 : Int nr = 0;
1241 527347 : for (Int irow=0; irow!=s.nRows(); irow++) {
1242 500650 : if (s.flagRow()(irow))
1243 382231 : continue;
1244 : Int iant;
1245 498997 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
1246 498997 : if (a1 == a2) {
1247 50223 : continue;
1248 : }
1249 448774 : else if (a1 == refant_) {
1250 118419 : iant = a2;
1251 : }
1252 330355 : else if (a2 == refant_) {
1253 0 : iant = a1;
1254 : }
1255 : else {
1256 330355 : continue;
1257 : }
1258 : // OK, we're not skipping this one so we have to do something.
1259 :
1260 : // v has shape (nelems, ?, nrows, nchannels)
1261 118419 : Cube<Complex> v = s.visCubeCorrected();
1262 118419 : Cube<Float> w = s.weightSpectrum();
1263 118419 : Cube<Bool> fl = s.flagCube();
1264 118419 : Int spw = s.spectralWindow()(0);
1265 118419 : Int f_index = gm_.swStartIndex(spw);
1266 118419 : Int t_index = gm_.getTimeIndex(s.time()(0));
1267 118419 : Int spwchans = gm_.getNChans(spw);
1268 118419 : IPosition start(4, 0, iant, t_index, f_index);
1269 118419 : IPosition stop(4, nCorr_, 1, 1, spwchans);
1270 118419 : IPosition stride(4, 1, 1, 1, 1);
1271 118419 : Slicer sl1(start, stop, stride, Slicer::endIsLength);
1272 236838 : Slicer sl2(IPosition(3, 0, 0, irow),
1273 236838 : IPosition(3, nCorr_, spwchans, 1),
1274 236838 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
1275 :
1276 236838 : Slicer flagSlice(IPosition(3, 0, 0, irow),
1277 236838 : IPosition(3, nCorr_, spwchans, 1),
1278 236838 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
1279 118419 : nr++;
1280 : if (DEVDEBUG && false) {
1281 : cerr << "nr " << nr
1282 : << " irow " << endl
1283 : << "Vpad shape " << Vpad_.shape() << endl
1284 : << "v shape " << v.shape() << endl
1285 : << "sl2 " << sl2 << endl
1286 : << "sl1 " << sl1 << endl
1287 : << "flagSlice " << flagSlice << endl;
1288 : }
1289 118419 : Array<Complex> rhs = v(sl2).nonDegenerate(1);
1290 118419 : Array<Float> weights = w(sl2).nonDegenerate(1);
1291 :
1292 118419 : unitize(rhs);
1293 118419 : Vpad_(sl1).nonDegenerate(1) = rhs * weights;
1294 :
1295 118419 : Array<Bool> flagged(fl(flagSlice).nonDegenerate(1));
1296 : // Zero flagged entries.
1297 118419 : Vpad_(sl1).nonDegenerate(1)(flagged) = Complex(0.0);
1298 :
1299 118419 : if (!allTrue(flagged)) {
1300 301639 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
1301 183220 : IPosition p(2, icorr, iant);
1302 183220 : Bool actually = false;
1303 183220 : activeAntennas_[icorr].insert(iant);
1304 3299236 : for (Int ichan=0; ichan != (Int) spwchans; ichan++) {
1305 3116016 : IPosition pchan(2, icorr, ichan);
1306 3116016 : if (!flagged(pchan)) {
1307 3075068 : Float wv = weights(pchan);
1308 3075068 : if (wv < 0) {
1309 0 : cerr << "spwchans " << spwchans << endl;
1310 0 : cerr << "Negative weight << (" << wv << ") on row "
1311 0 : << irow << " baseline (" << a1 << ", " << a2 << ") "
1312 0 : << " channel " << ichan << endl;
1313 0 : cerr << "pchan " << pchan << endl;
1314 0 : cerr << "Weights " << weights << endl;
1315 : }
1316 3075068 : xcount_(p)++;
1317 3075068 : sumw_(p) += wv;
1318 3075068 : sumww_(p) += wv*wv;
1319 3075068 : actually = true;
1320 : }
1321 3116016 : }
1322 183220 : if (actually) {
1323 183220 : activeAntennas_[icorr].insert(iant);
1324 183220 : goodRows++;
1325 : }
1326 183220 : }
1327 : }
1328 : if (DEVDEBUG && 0) {
1329 : cerr << "flagged " << flagged << endl;
1330 : cerr << "flagSlice " << flagSlice << endl
1331 : << "fl.shape() " << fl.shape() << endl
1332 : << "Vpad_.shape() " << Vpad_.shape() << endl
1333 : << "flagged.shape() " << flagged.shape() << endl
1334 : << "sl1 " << sl1 << endl;
1335 : }
1336 118419 : }
1337 : }
1338 : if (DEVDEBUG) {
1339 : cerr << "In DelayRateFFTConcat::DelayRateFFTConcat " << endl;
1340 : printActive();
1341 : cerr << "sumw_ " << sumw_ << endl;
1342 : cerr << "Constructed a DelayRateFFTConcat object." << endl;
1343 : cerr << "totalRows " << totalRows << endl;
1344 : cerr << "goodRows " << goodRows << endl;
1345 : }
1346 :
1347 217 : }
1348 :
1349 0 : DelayRateFFTConcat::DelayRateFFTConcat(Array<Complex>& data, Int nPadFactor, Float f0, Float df, Float dt, SDBList& s,
1350 0 : Array<Double>& delayWindow, Array<Double>& rateWindow) :
1351 : DelayRateFFT(s, 0, delayWindow, rateWindow),
1352 0 : gm_(s),
1353 0 : nPadFactor_(nPadFactor),
1354 0 : Vpad_() {
1355 0 : dt_ = dt;
1356 0 : f0_ = f0;
1357 0 : df_ = df;
1358 0 : IPosition shape = data.shape();
1359 0 : nCorr_ = shape(0);
1360 0 : nElem_ = shape(1);
1361 0 : nt_ = shape(2);
1362 0 : nChan_ = shape(3);
1363 0 : nPadT_ = nPadFactor_*nt_;
1364 0 : nPadChan_ = nPadFactor_*nChan_;
1365 0 : IPosition paddedDataSize(4, nCorr_, nElem_, nPadT_, nPadChan_);
1366 0 : Vpad_.resize(paddedDataSize);
1367 :
1368 0 : IPosition start(4, 0, 0, 0, 0);
1369 0 : IPosition stop(4, nCorr_, nElem_, nt_, nChan_);
1370 0 : IPosition stride(4, 1, 1, 1, 1);
1371 0 : Slicer sl1(start, stop, stride, Slicer::endIsLength);
1372 0 : Vpad_(sl1) = data;
1373 :
1374 0 : unitize(Vpad_);
1375 :
1376 0 : }
1377 :
1378 : void
1379 217 : DelayRateFFTConcat::FFT() {
1380 : if (DEVDEBUG) {
1381 : cerr << "DelayRateFFTConcat::FFT()" << endl;
1382 : }
1383 : // Axes are 0: correlation (i.e., hand of polarization), 1: antenna, 2: time, 3: channel
1384 217 : Vector<Bool> ax(4, false);
1385 217 : ax(2) = true;
1386 217 : ax(3) = true;
1387 : // Also copied from DelayFFT in KJones.
1388 : // we make a copy to FFT in place.
1389 : if (DEVDEBUG) {
1390 : cerr << "Vpad_.shape() " << Vpad_.shape() << endl;
1391 : }
1392 217 : ArrayLattice<Complex> c(Vpad_);
1393 217 : LatticeFFT::cfft0(c, ax, true);
1394 : if (DEVDEBUG) {
1395 : cerr << "FFT transformed" << endl;
1396 : }
1397 217 : }
1398 :
1399 : std::pair<Bool, Float>
1400 4076 : DelayRateFFTConcat::xinterp(Float alo, Float amax, Float ahi) {
1401 4076 : if (amax > 0.0 && alo == amax && amax == ahi)
1402 0 : return std::make_pair(true, 0.5);
1403 :
1404 4076 : Float denom(alo-2.0*amax+ahi);
1405 4076 : Bool cond = amax>0.0 && abs(denom)>0.0 ;
1406 4076 : Float fpk = cond ? 0.5-(ahi-amax)/denom : 0.0;
1407 4076 : return std::make_pair(cond, fpk);
1408 : }
1409 :
1410 : void
1411 217 : DelayRateFFTConcat::searchPeak() {
1412 : if (DEVDEBUG) {
1413 : cerr << "DelayRateFFTConcat::searchPeak()" << endl;
1414 : }
1415 :
1416 : // Recall param_ -> [phase, delay, rate] for each correlation
1417 217 : param_.resize(3*nCorr_, nElem_); // This might be better done elsewhere.
1418 217 : param_.set(0.0);
1419 217 : flag_.resize(3*nCorr_, nElem_);
1420 217 : flag_.set(true); // all flagged initially
1421 : if (DEVDEBUG) {
1422 : cerr << "nt_ " << nt_ << " nPadChan_ " << nPadChan_ << endl;
1423 : cerr << "Vpad_.shape() " << Vpad_.shape() << endl;
1424 : cerr << "delayWindow_ " << delayWindow_ << endl;
1425 :
1426 : }
1427 :
1428 567 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
1429 350 : flag_(icorr*3 + 0, refant()) = false;
1430 350 : flag_(icorr*3 + 1, refant()) = false;
1431 350 : flag_(icorr*3 + 2, refant()) = false;
1432 2738 : for (Int ielem=0; ielem<nElem_; ++ielem) {
1433 2388 : if (ielem==refant()) {
1434 350 : continue;
1435 : }
1436 : // NB: Time, Channel
1437 : // And once again we fail at slicing
1438 2038 : IPosition start(4, icorr, ielem, 0, 0);
1439 2038 : IPosition stop(4, 1, 1, nPadT_, nPadChan_);
1440 2038 : IPosition step(4, 1, 1, 1, 1);
1441 2038 : Slicer sl(start, stop, step, Slicer::endIsLength);
1442 4076 : Matrix<Complex> aS = Vpad_(sl).nonDegenerate();
1443 2038 : Int sgn = (ielem < refant()) ? 1 : -1;
1444 :
1445 : // Below is the gory details for turning delay window into index range
1446 2038 : Double bw = Float(nPadChan_)*df_;
1447 2038 : Double d0 = sgn*delayWindow_(IPosition(1, 0));
1448 2038 : Double d1 = sgn*delayWindow_(IPosition(1, 1));
1449 2038 : if (d0 > d1) std::swap(d0, d1);
1450 2038 : d0 = max(d0, -0.5/df_);
1451 2038 : d1 = min(d1, (0.5-1/nPadChan_)/df_);
1452 :
1453 : // It's simpler to keep the ranges as signed integers and
1454 : // handle the wrapping of the FFT in the loop over
1455 : // indices. Recall that the FFT result returned has indices
1456 : // that run from 0 to nPadChan_/2 -1 and then from
1457 : // -nPadChan/2 to -1, so far as our delay is concerned.
1458 2038 : Int i0 = bw*d0;
1459 2038 : Int i1 = bw*d1;
1460 2038 : if (i1==i0) i1++;
1461 : // Now for the gory details for turning rate window into index range
1462 2038 : Double width = nPadT_*dt_*1e9*f0_;
1463 2038 : Double r0 = sgn*rateWindow_(IPosition(1,0));
1464 2038 : Double r1 = sgn*rateWindow_(IPosition(1,1));
1465 2038 : if (r0 > r1) std::swap(r0, r1);
1466 2038 : r0 = max(r0, -0.5/(dt_*1e9*f0_));
1467 2038 : r1 = min(r1, (0.5 - 1/nPadT_)/(dt_*1e9*f0_));
1468 :
1469 2038 : Int j0 = width*r0;
1470 2038 : Int j1 = width*r1;
1471 2038 : if (j1==j0) j1++;
1472 : if (DEVDEBUG) {
1473 : cerr << "Checking the windows for delay and rate search." << endl;
1474 : cerr << "bw " << bw << endl;
1475 : cerr << "d0 " << d0 << " d1 " << d1 << endl;
1476 : cerr << "i0 " << i0 << " i1 " << i1 << endl;
1477 : cerr << "r0 " << r0 << " r1 " << r1 << endl;
1478 : cerr << "j0 " << j0 << " j1 " << j1 << endl;
1479 : }
1480 2038 : Matrix<Float> amp(amplitude(aS));
1481 2038 : Int ipkch(0);
1482 2038 : Int ipkt(0);
1483 2038 : Float amax(-1.0);
1484 : // Unlike KJones we have to iterate in time too
1485 1645356 : for (Int itime0=j0; itime0 != j1; itime0++) {
1486 1643318 : Int itime = (itime0 < 0) ? itime0 + nPadT_ : itime0;
1487 223883766 : for (Int ich0=i0; ich0 != i1; ich0++) {
1488 222240448 : Int ich = (ich0 < 0) ? ich0 + nPadChan_ : ich0;
1489 : // cerr << "Gridpoint " << itime << ", " << ich << "->" << amp(itime, ich) << endl;
1490 222240448 : if (amp(itime, ich) > amax) {
1491 196045 : ipkch = ich;
1492 196045 : ipkt = itime;
1493 196045 : amax=amp(itime, ich);
1494 : }
1495 : }
1496 : }
1497 : // Finished grovelling. Now we have the location of the
1498 : // maximum amplitude.
1499 2038 : Float alo_ch = amp(ipkt, (ipkch > 0) ? ipkch-1 : nPadChan_-1);
1500 2038 : Float ahi_ch = amp(ipkt, ipkch<(nPadChan_-1) ? ipkch+1 : 0);
1501 2038 : std::pair<Bool, Float> maybeFpkch = xinterp(alo_ch, amax, ahi_ch);
1502 : // We handle wrapping while looking for neighbours
1503 2038 : Float alo_t = amp(ipkt > 0 ? ipkt-1 : nPadT_ -1, ipkch);
1504 2038 : Float ahi_t = amp(ipkt < (nPadT_ -1) ? ipkt+1 : 0, ipkch);
1505 : if (DEVDEBUG) {
1506 : cerr << "For element " << ielem << endl;
1507 : cerr << "In channel dimension ipkch " << ipkch << " alo " << alo_ch
1508 : << " amax " << amax << " ahi " << ahi_ch << endl;
1509 : cerr << "In time dimension ipkt " << ipkt << " alo " << alo_t
1510 : << " amax " << amax << " ahi " << ahi_t << endl;
1511 : }
1512 2038 : std::pair<Bool, Float> maybeFpkt = xinterp(alo_t, amax, ahi_t);
1513 :
1514 2038 : if (maybeFpkch.first and maybeFpkt.first) {
1515 : // Phase
1516 1801 : Complex c = aS(ipkt, ipkch);
1517 1801 : Float phase = arg(c);
1518 1801 : param_(icorr*3 + 0, ielem) = sgn*phase;
1519 1801 : Float delay = (ipkch)/Float(nPadChan_);
1520 1801 : if (delay > 0.5) delay -= 1.0; // fold
1521 1801 : delay /= df_; // nsec
1522 1801 : param_(icorr*3 + 1, ielem) = sgn*delay; //
1523 1801 : Double rate = (ipkt)/Float(nPadT_);
1524 1801 : if (rate > 0.5) rate -= 1.0;
1525 1801 : Double rate0 = rate/dt_;
1526 1801 : Double rate1 = rate0/(1e9 * f0_);
1527 :
1528 1801 : param_(icorr*3 + 2, ielem) = Float(sgn*rate1);
1529 : if (DEVDEBUG) {
1530 : cerr << "maybeFpkch.second=" << maybeFpkch.second
1531 : << ", df_= " << df_
1532 : << " fpkch=" << (ipkch + maybeFpkch.second) << endl;
1533 : cerr << " maybeFpkt.second=" << maybeFpkt.second
1534 : << " rate0=" << rate
1535 : << " 1e9 * f0_=" << 1e9 * f0_
1536 : << ", dt_=" << dt_
1537 : << " fpkt=" << (ipkt + maybeFpkt.second) << endl;
1538 :
1539 : }
1540 : if (DEVDEBUG) {
1541 : cerr << "Found peak for element " << ielem << " correlation " << icorr
1542 : << " ipkt=" << ipkt << "/" << nPadT_ << ", ipkch=" << ipkch << "/" << nPadChan_
1543 : << " peak=" << amax
1544 : << "; delay " << delay << ", rate " << rate
1545 : << ", phase " << arg(c) << " sign= " << sgn << endl;
1546 : }
1547 : // Set 3 flags.
1548 1801 : flag_(icorr*3 + 0, ielem)=false;
1549 1801 : flag_(icorr*3 + 1, ielem)=false;
1550 1801 : flag_(icorr*3 + 2, ielem)=false;
1551 1801 : }
1552 : else {
1553 : if (DEVDEBUG) {
1554 : cerr << "No peak in 2D FFT for element " << ielem << " correlation " << icorr << endl;
1555 : }
1556 : }
1557 2038 : }
1558 : }
1559 217 : }
1560 :
1561 :
1562 : Float
1563 1802 : DelayRateFFTConcat::snr(Int icorr, Int ielem, Float delay, Float rate) {
1564 : if (DEVDEBUG) {
1565 : cerr << "DelayRateFFTConcat::snr"<< endl;
1566 : }
1567 : // We calculate a signal-to-noise ration for the 2D FFT fringefit
1568 : // using a formula transcribed from AIPS FRING.
1569 : //
1570 : // Have to convert delay and rate back into indices on the padded 2D grid.
1571 1802 : Int sgn = (ielem < refant()) ? 1 : -1;
1572 1802 : delay *= sgn*df_;
1573 1802 : if (delay < 0.0) delay += 1;
1574 1802 : Int ichan = Int(delay*nPadChan_ + 0.5);
1575 1802 : if (ichan == nPadChan_) ichan = 0;
1576 :
1577 1802 : rate *= sgn*1e9 * f0_;
1578 1802 : rate *= dt_;
1579 1802 : if (rate < 0.0) rate += 1;
1580 1802 : Int itime = Int(rate*nPadT_ + 0.5);
1581 1802 : if (itime == nPadT_) itime = 0;
1582 : // What about flags? If the datapoint closest to the computed
1583 : // delay and rate values is flagged we probably shouldn't use
1584 : // it, but what *should* we use?
1585 1802 : IPosition ipos(4, icorr, ielem, itime, ichan);
1586 1802 : IPosition p(2, icorr, ielem);
1587 1802 : Complex v = Vpad_(ipos);
1588 1802 : Float peak = abs(v);
1589 1802 : if (peak > 0.999*sumw_(p)) peak=0.999*sumw_(p);
1590 : // xcount is number of data points for baseline to ielem
1591 : // sumw is sum of weights,
1592 : // sumww is sum of squares of weights
1593 :
1594 : Float cwt;
1595 1802 : if (fabs(sumw_(p))<FLT_EPSILON) {
1596 1 : cwt = 0;
1597 : if (DEVDEBUG) {
1598 : cerr << "Correlation " << icorr << " antenna " << ielem << ": sum of weights is zero." << endl;
1599 : }
1600 : } else {
1601 1801 : Float x = M_PI/2*peak/sumw_(p);
1602 : // The magic numbers in the following formula are from AIPS FRING
1603 1801 : cwt = (pow(tan(x), 1.163) * sqrt(sumw_(p)/sqrt(sumww_(p)/xcount_(p))));
1604 : if (DEVDEBUG) {
1605 : cerr << "Correlation " << icorr << " antenna " << ielem << " ipos " << ipos
1606 : << " peak=" << peak << "; xang=" << x << "; xcount=" << xcount_(p) << "; sumw=" << sumw_(p) << "; sumww=" << sumww_(p)
1607 : << " snr " << cwt << endl;
1608 : }
1609 : }
1610 1802 : return cwt;
1611 1802 : }
1612 :
1613 :
1614 : } // end namespace
|
