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 112487 : static void unitize(Array<Complex>& vC)
76 : {
77 112487 : Array<Float> vCa(amplitude(vC));
78 : // Divide by non-zero amps
79 112487 : vCa(vCa<FLT_EPSILON)=1.0;
80 112487 : vC /= vCa;
81 112487 : }
82 :
83 0 : SDBListGridManagerCombo::SDBListGridManagerCombo(SDBList& sdbs) :
84 : SDBListGridManager(sdbs),
85 0 : nchan_(0)
86 : {
87 0 : std::set<Double> fmaxes_;
88 0 : std::set<Double> fmins_;
89 :
90 0 : 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 0 : for (Int i=0; i != sdbs_.nSDB(); i++) {
100 0 : SolveDataBuffer& sdb = sdbs_(i);
101 0 : Int pspw = sdb.spectralWindow()(0);
102 0 : Double t = sdbs_(i).time()(0);
103 0 : times_.insert(t);
104 0 : if (spwins_.find(pspw) == spwins_.end()) {
105 0 : spwins_.insert(pspw);
106 0 : const Vector<Double>& fs = sdb.freqs();
107 0 : df_ = fs[1] - fs[0];
108 0 : pspwIdToFreqMap_[pspw] = &(sdb.freqs());
109 0 : nchan_ = max(nchan_, sdb.nChannels());
110 0 : fmaxes_.insert(fs(nchan_-1));
111 0 : fmins_.insert(fs(0));
112 : } else {
113 0 : continue;
114 : }
115 : }
116 0 : nt_ = times_.size();
117 : // nt_ = sdbs_.nSDB()/spwins_.size();
118 0 : tmin_ = *(times_.begin());
119 0 : tmax_ = *(times_.rbegin());
120 0 : dt_ = (tmax_ - tmin_)/(nt_ - 1);
121 0 : if (nchan_ == 1) {
122 0 : df_ = 1;
123 0 : return;
124 : }
125 0 : size_t i = 0;
126 0 : for (auto p=spwins_.begin(); p!=spwins_.end(); p++) {
127 0 : spwPMap_[*p] = i;
128 0 : i++;
129 : }
130 0 : }
131 :
132 :
133 :
134 :
135 : Float
136 0 : SDBListGridManagerCombo::getRefFreqFromLSPW(Int lspw) {
137 0 : auto p = spwins_.begin();
138 0 : std::advance(p, lspw);
139 0 : Int pspw = *p;
140 0 : const Vector<Double>& fs = *(pspwIdToFreqMap_[pspw]);
141 0 : return fs[0];
142 : }
143 :
144 219 : DelayRateFFT::DelayRateFFT(SDBList& sdbs, Int refant, Array<Double>& delayWindow, Array<Double>& rateWindow) :
145 438 : f0_(sdbs.centroidFreq() / 1.e9), // GHz, for delayrate calc
146 219 : refant_(refant),
147 219 : delayWindow_(delayWindow),
148 219 : rateWindow_(rateWindow),
149 219 : Vall_(),
150 219 : xcount_(),
151 219 : sumw_(),
152 219 : sumww_(),
153 219 : activeAntennas_(),
154 657 : allActiveAntennas_()
155 219 : {}
156 :
157 : DelayRateFFT *
158 219 : DelayRateFFT::makeAChild(int concat, SDBList& sdbs,
159 : Int refant,
160 : Array<Double>& delayWindow_,
161 : Array<Double>& rateWindow_)
162 : {
163 219 : if (concat) {
164 219 : return new DelayRateFFTConcat(sdbs, refant, delayWindow_, rateWindow_);
165 : } else {
166 0 : 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 356 : void DelayRateFFT::removeAntennasCorrelation(Int icor, std::set< Int > s) {
190 356 : std::set< Int > & as = activeAntennas_.find(icor)->second;
191 356 : for (std::set< Int >::iterator it=s.begin(); it!=s.end(); it++) {
192 0 : as.erase(*it);
193 : }
194 356 : }
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 0 : DelayRateFFTCombo::DelayRateFFTCombo(SDBList& sdbs, Int refant, Array<Double>& delayWindow,
213 0 : Array<Double>& rateWindow) :
214 : DelayRateFFT(sdbs, refant, delayWindow, rateWindow),
215 0 : gm_(sdbs)
216 : {
217 : if (DEVDEBUG) {
218 : gm_.describe();
219 : }
220 0 : nt_ = gm_.nt_;
221 0 : nChan_ = gm_.nchan_;
222 :
223 : // We might want to pad these too
224 0 : nPadFactor_ = 4;
225 0 : nPadT_ = nPadFactor_*nt_;
226 0 : nPadChan_ = nPadFactor_*nChan_;
227 0 : nspw_ = gm_.nSPW();
228 0 : dt_ = gm_.dt_;
229 0 : df_ = gm_.df_ / 1.e9;
230 :
231 0 : 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 0 : IPosition ds = delayWindow_.shape();
235 0 : if (ds.size()!=1 || ds.nelements() != 1) {
236 0 : throw AipsError("delaywindow must be a list of length 2.");
237 : }
238 0 : IPosition rs = rateWindow_.shape();
239 0 : if (rs.size()!=1 || rs.nelements() != 1) {
240 0 : throw AipsError("ratewindow must be a list of length 2.");
241 : }
242 0 : Int nCorrOrig(sdbs(0).nCorrelations());
243 0 : nCorr_ = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
244 0 : for (Int i=0; i<nCorr_; i++) {
245 0 : 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 0 : Int corrStep = (nCorrOrig > 2 ? 3 : 1); // step for p-hands
251 0 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
252 0 : SolveDataBuffer& s(sdbs(ibuf));
253 0 : for (Int irow=0; irow!=s.nRows(); irow++) {
254 0 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
255 0 : allActiveAntennas_.insert(a1);
256 0 : allActiveAntennas_.insert(a2);
257 : }
258 : }
259 0 : nElem_ = 1 + *(allActiveAntennas_.rbegin()) ;
260 :
261 0 : IPosition aggregateDim(2, nCorr_, nElem_);
262 0 : xcount_.resize(aggregateDim);
263 0 : sumw_.resize(aggregateDim);
264 0 : sumww_.resize(aggregateDim);
265 0 : peak_.resize(aggregateDim);
266 :
267 0 : xcount_ = 0;
268 0 : sumw_ = 0.0;
269 0 : sumww_ = 0.0;
270 0 : IPosition dataSize(5, nCorr_, nElem_, nspw_, nPadT_, nPadChan_);
271 0 : Vall_.resize(dataSize);
272 0 : Int totalRows = 0;
273 0 : Int goodRows = 0;
274 :
275 : if (DEVDEBUG) {
276 : cerr << "DelayRateFFTCombo constructor: Here" << endl;
277 : }
278 0 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
279 0 : SolveDataBuffer& s(sdbs(ibuf));
280 0 : totalRows += s.nRows();
281 0 : if (!s.Ok())
282 0 : continue;
283 :
284 0 : Int nr = 0;
285 0 : for (Int irow=0; irow!=s.nRows(); irow++) {
286 0 : if (s.flagRow()(irow))
287 0 : continue;
288 : Int iant;
289 0 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
290 0 : if (a1 == a2) { continue; } // we don't do autocorrelations
291 0 : else if (a1 == refant_) { iant = a2; }
292 0 : else if (a2 == refant_) { iant = a1; }
293 0 : else { continue; } // not a baseline to reference antenna
294 : // v has shape (nelems, ?, nrows, nchannels); can't be a reference.
295 0 : Cube<Complex> v = s.visCubeCorrected();
296 0 : const Cube<Float>& w( s.weightSpectrum() );
297 0 : const Cube<Bool>& fl( s.flagCube() );
298 0 : Int pspw = s.spectralWindow()(0);
299 0 : Int ispw = gm_.getLSPW(pspw);
300 0 : 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 0 : IPosition start (5, 0, iant, ispw, t_index, 0);
308 0 : IPosition stop (5, nCorr_, 1, 1, 1, nChan_);
309 0 : IPosition stride(5, 1, 1, 1, 1, 1);
310 0 : 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 0 : Slicer source_slice(IPosition(3, 0, 0, irow),
328 0 : IPosition(3, nCorr_, nChan_, 1),
329 0 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
330 :
331 0 : Slicer flagSlice(IPosition(3, 0, 0, irow),
332 0 : IPosition(3, nCorr_, nChan_, 1),
333 0 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
334 0 : nr++;
335 0 : Array<Complex>rhs( v(source_slice).nonDegenerate(1) );
336 0 : const Array<Float>& weights( w(source_slice).nonDegenerate(1) );
337 0 : unitize(rhs);
338 0 : Vall_(target_slice).nonDegenerate(1) = rhs * weights;
339 0 : const Array<Bool>& flagged( fl(flagSlice).nonDegenerate(1) );
340 :
341 : // Zero flagged entries.
342 0 : Vall_(target_slice).nonDegenerate(1)(flagged) = Complex(0.0);
343 :
344 0 : if (!allTrue(flagged)) {
345 0 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
346 : // Book keeping now done per spw!
347 0 : IPosition p(3, icorr, iant, ispw);
348 0 : Bool actually = false;
349 0 : activeAntennas_[icorr].insert(iant);
350 0 : for (Int ichan=0; ichan != (Int) nChan_; ichan++) {
351 0 : IPosition pchan(2, icorr, ichan);
352 0 : if (!flagged(pchan)) {
353 0 : Float wv = weights(pchan);
354 0 : xcount_(p)++;
355 0 : sumw_(p) += wv;
356 0 : sumww_(p) += wv*wv;
357 0 : actually = true;
358 : }
359 0 : }
360 0 : if (actually) {
361 0 : activeAntennas_[icorr].insert(iant);
362 0 : goodRows++;
363 : }
364 0 : }
365 : }
366 0 : }
367 : }
368 0 : }
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 0 : 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 0 : Vector<Bool> ax(5, false);
406 0 : ax(3) = true;
407 0 : ax(4) = true;
408 : // Also copied from DelayFFT in KJones.
409 : // we make a copy to FFT in place.
410 0 : ArrayLattice<Complex> c(Vall_);
411 0 : LatticeFFT::cfft0(c, ax, true);
412 : if (DEVDEBUG) {
413 : cerr << "FFT transformed" << endl;
414 : }
415 0 : }
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 0 : DelayRateFFTCombo::searchPeak() {
423 : if (DEVDEBUG) {
424 : cerr << "DelayRateFFTCombo::searchPeak()" << endl;
425 : }
426 :
427 : // Recall param_ -> [phase, delay, rate] for each correlation
428 0 : param_.resize(3*nCorr_, nElem_); // This might be better done elsewhere.
429 0 : param_.set(0.0);
430 0 : flag_.resize(3*nCorr_, nElem_);
431 0 : flag_.set(true); // all flagged initially
432 :
433 0 : Double bw = Float(nChan_)*df_;
434 0 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
435 0 : flag_(icorr*3 + 0, refant()) = false;
436 0 : flag_(icorr*3 + 1, refant()) = false;
437 0 : flag_(icorr*3 + 2, refant()) = false;
438 0 : for (Int ielem=0; ielem<nElem_; ++ielem) {
439 0 : if (ielem==refant()) {
440 0 : continue;
441 : }
442 0 : if (activeAntennas_[icorr].find(ielem)==activeAntennas_[icorr].end()) {
443 0 : continue;
444 : }
445 :
446 : // Below is the gory details for turning delay window into index range
447 0 : Int sgn = (ielem < refant()) ? 1 : -1;
448 0 : Double d0 = sgn*delayWindow_(IPosition(1, 0));
449 0 : Double d1 = sgn*delayWindow_(IPosition(1, 1));
450 0 : if (d0 > d1) std::swap(d0, d1);
451 0 : d0 = max(d0, -0.5/df_);
452 0 : 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 0 : Int i0 = bw*d0;
460 0 : Int i1 = bw*d1;
461 0 : i0 *= nPadFactor_;
462 0 : i1 *= nPadFactor_;
463 0 : if (i1==i0) i1++;
464 : // Now for the gory details for turning rate window into index range
465 0 : Double width = nt_*dt_*1e9*f0_;
466 0 : Double r0 = sgn*rateWindow_(IPosition(1,0));
467 0 : Double r1 = sgn*rateWindow_(IPosition(1,1));
468 0 : if (r0 > r1) std::swap(r0, r1);
469 0 : r0 = max(r0, -0.5/(dt_*1e9*f0_));
470 0 : r1 = min(r1, (0.5 - 1/nt_)/(dt_*1e9*f0_));
471 :
472 0 : Int j0 = width*r0;
473 0 : Int j1 = width*r1;
474 0 : j0 *= nPadFactor_;
475 0 : j1 *= nPadFactor_;
476 0 : if (j1==j0) j1++;
477 : Array<Complex> blVis(
478 0 : Vall_(Slicer(
479 0 : IPosition(5, icorr, ielem, 0, 0, 0),
480 0 : IPosition(5, 1, 1, nspw_, nPadT_, nPadChan_),
481 0 : IPosition(5, 1, 1, 1, 1, 1),
482 0 : Slicer::endIsLength)).nonDegenerate(IPosition(1,2)));
483 :
484 0 : Vector<Float> offsets(nspw_);
485 0 : Double bw = gm_.nchan_ * gm_.df_;
486 0 : for (size_t lspw=0; lspw!=size_t(nspw_); lspw++) {
487 0 : Double f = gm_.getRefFreqFromLSPW(lspw);
488 0 : Double f0 = gm_.getRefFreqFromLSPW(0);
489 0 : offsets(lspw) = (f - f0)/bw;
490 : }
491 :
492 : // Get the DFT estimates of peak location instead
493 0 : tuple<Double, Double, Double, Double> mp = multibandFFTs(blVis , offsets);
494 0 : Double mpkt = std::get<0>(mp);
495 0 : Double mpkch = std::get<1>(mp);
496 0 : Double mpeak = std::get<2>(mp);
497 0 : Double marg = std::get<3>(mp);
498 :
499 0 : Complex c0 = dotWithOffsets(blVis, offsets, mpkt, mpkch);
500 : if (DEVDEBUG) {
501 : cerr << "DelayRateFFTCombo: abs(c1) " << abs(c0) << endl;
502 : }
503 0 : tuple<Double, Double, Double, Double> p = refineSearch(blVis, offsets, mpkt, mpkch);
504 0 : Double pkt = std::get<0>(p);
505 0 : Double pkch = std::get<1>(p);
506 0 : Double peak = std::get<2>(p);
507 0 : 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 0 : peak_(IPosition(2, icorr, ielem)) = peak;
517 0 : param_(icorr*3 + 0, ielem) = sgn*phase;
518 0 : Double delay = (pkch)/Double(nChan_);
519 0 : if (delay > 0.5) delay -= 1.0; // fold
520 0 : delay /= df_; // nsec
521 0 : param_(icorr*3 + 1, ielem) = sgn*delay;
522 0 : Double rate = (pkt)/Double(nt_);
523 0 : if (rate > 0.5) rate -= 1.0;
524 0 : Double rate0 = rate/dt_;
525 0 : Double rate1 = rate0/(1e9 * f0_);
526 0 : 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 0 : flag_(icorr*3 + 0, ielem)=false;
532 0 : flag_(icorr*3 + 1, ielem)=false;
533 0 : 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 0 : }
539 : }
540 0 : }
541 :
542 :
543 : tuple<Double, Double, Double, Double>
544 0 : 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 0 : IPosition shp = ffts.shape();
559 0 : Int nspw = shp(0);
560 0 : Int nt = shp(1);
561 0 : Int nchan = shp(2);
562 :
563 0 : Int binScale = 1;
564 0 : size_t nbins= size_t(binScale*max(offsets)+1+0.5);
565 0 : Vector<Double> freqs(nbins);
566 0 : freqs = 0;
567 0 : for (size_t k=0; k!=nbins; k++) {
568 0 : 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 0 : Array<Complex> X(IPosition(3, nbins, nt, nchan));
573 0 : X = 0;
574 :
575 0 : for (size_t ispw=0; ispw!=size_t(nspw); ispw++) {
576 0 : for (size_t it=0; it != nt; it++) {
577 0 : for (size_t ichan=0; ichan != nchan; ichan++) {
578 0 : Double x = abs(ffts(IPosition(3, ispw, 0, 0)));
579 0 : if (std::isnan(x)) {cerr << "Skipping FFT " << ispw << endl;
580 0 : continue;
581 : }
582 0 : for (size_t ibin=0; ibin!=nbins; ibin++) {
583 0 : Complex rotation = exp(-Complex(0,1)*Complex(C::_2pi*offsets(ispw)*freqs(ibin)));
584 0 : 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 0 : Int itmax = -1;
633 0 : Int ichanmax = -1;
634 0 : Int ibinmax = -1;
635 0 : Double zmax = -1.0;
636 0 : for (Int ibin=0; ibin != nbins; ibin++) {
637 0 : for (Int it = 0; it != nt; it++) {
638 0 : for (Int ichan=0; ichan != nchan; ichan++) {
639 0 : Complex c = X(IPosition(3, ibin, it, ichan));
640 0 : Double a = abs(c);
641 0 : if (a>zmax) {
642 0 : itmax = it;
643 0 : ichanmax = ichan;
644 0 : ibinmax = ibin;
645 0 : 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 0 : Complex c = X(IPosition(3, ibinmax, itmax, ichanmax));
659 0 : Double dpkch = freqs(ibinmax);
660 :
661 0 : tuple<Double, Double, Double, Double> t;
662 0 : 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 0 : return t;
671 0 : }
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(C::_2pi*offsets(n)*freqs(k)));
724 : // X[k] += phasors[n]*exp(-Complex(0,1)*Complex(C::_2pi*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 0 : dotWithOffsets(const Cube<Complex>& ft, const Vector<Float>& offsets, double k, double l) {
751 0 : Int nspw = ft.nrow();
752 0 : Int ni = ft.ncolumn();
753 0 : Int nj = ft.nplane();
754 :
755 0 : if (k<0) k += (ni);
756 0 : if (l<0) l += (nj);
757 0 : Complex p(0.0, 0.0);
758 0 : for (Int s=0; s!=nspw; s++) {
759 0 : const Matrix<Complex>& ft_s = ft.yzPlane(s);
760 0 : Double f_off = offsets(s);
761 0 : Complex r = dotMatrixWithModel(ft_s, k, l, f_off);
762 0 : p += r;
763 0 : }
764 0 : 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 0 : c_peak_fn(const gsl_vector *x, void *vparams) {
772 0 : std::pair< Cube<Complex> const *, Vector<Float> const * > *params =
773 : (std::pair< Cube<Complex> const *, Vector<Float> const * > *)(vparams );
774 :
775 0 : double k = gsl_vector_get(x, 0);
776 0 : double l = gsl_vector_get(x, 1);
777 :
778 0 : const Cube<Complex>& ft ( *(params->first ));
779 0 : const Vector<Float>& offsets ( *(params->second));
780 0 : 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 0 : my_peak_fn(const gsl_vector *x, void *vparams) {
788 0 : Complex p = c_peak_fn(x, vparams);
789 0 : 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 0 : 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 0 : const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
825 : /* Starting point */
826 0 : gsl_vector *x = gsl_vector_alloc (2);
827 0 : gsl_vector_set(x, 0, pkt);
828 0 : 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 0 : gsl_vector* steps = gsl_vector_alloc (2);
835 0 : gsl_vector_set_all(steps, 0.01);
836 :
837 : /* Initialize method and iterate */
838 :
839 : gsl_multimin_function minex_func;
840 0 : std::pair< Cube<Complex> const * , Vector<Float> const * > par = make_pair(&ft, &offsets);
841 :
842 0 : minex_func.n = 2;
843 0 : minex_func.f = my_peak_fn;
844 0 : minex_func.params = ∥
845 :
846 0 : gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, 2);
847 0 : gsl_multimin_fminimizer_set (s, &minex_func, x, steps);
848 :
849 0 : float minstep = 1e-5;
850 : int status;
851 0 : size_t iter = 0;
852 : do {
853 : if (DEVDEBUG) {
854 : cerr << "Starting iteration" << endl;
855 : }
856 :
857 0 : iter++;
858 0 : status = gsl_multimin_fminimizer_iterate(s);
859 0 : if (status) break;
860 0 : double size = gsl_multimin_fminimizer_size(s);
861 0 : status = gsl_multimin_test_size(size, minstep);
862 0 : 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 0 : } while (status == GSL_CONTINUE && iter < 5);
874 : // while (status == GSL_CONTINUE && iter < 100);
875 0 : tuple<Double, Double, Double, Double> p;
876 : // if (status == GSL_SUCCESS) {
877 : if (1) {
878 : // Double peak = gsl_multimin_fminimizer_minimum(s);
879 0 : Double ipkt = gsl_vector_get(s->x, 0);
880 0 : Double ipkch = gsl_vector_get(s->x, 1);
881 0 : Complex c = c_peak_fn(s->x, &par);
882 0 : 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 0 : gsl_vector_free(steps);
887 0 : gsl_vector_free(x);
888 0 : gsl_multimin_fminimizer_free(s);
889 :
890 0 : return p;
891 : }
892 :
893 : Float
894 0 : 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 0 : IPosition p(2, icorr, ielem);
901 0 : Float peak = peak_(p);
902 0 : if (peak > 0.999*sumw_(p)) {
903 : //cerr << "Clipping peak for element " << ielem << " from " << peak << " to " << 0.999*sumw_(p) << endl;
904 0 : 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 0 : if (fabs(sumw_(p))<FLT_EPSILON) {
911 0 : cwt = 0;
912 : if (DEVDEBUG) {
913 : cerr << "Correlation " << icorr << " antenna " << ielem << ": sum of weights is zero." << endl;
914 : }
915 : } else {
916 0 : Float x = C::pi/2*peak/sumw_(p);
917 : // The magic numbers in the following formula are from AIPS FRING
918 0 : 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 0 : return cwt;
926 0 : }
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 = C::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 0 : dotMatrixWithModel(const Matrix<Complex>& data, Double k, Double l, Double offset)
993 : {
994 0 : size_t ni = data.nrow();
995 0 : size_t nj = data.ncolumn();
996 0 : 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 0 : Double eps = 1e-8;
1001 0 : Int k_int = floor(k);
1002 0 : Double k_del = k - k_int;
1003 0 : Bool k_flag = (fabs(k_del) < eps);
1004 0 : Int l_int = floor(l);
1005 0 : Double l_del = l - l_int;
1006 0 : Bool l_flag = (fabs(l_del) < eps);
1007 0 : Complex t0;
1008 0 : Complex t1;
1009 : // FIXME! We're messing with reversing signs
1010 0 : for (size_t i=0; i!=ni; i++) {
1011 0 : if (k_flag) {
1012 0 : t0 = Complex(i==k);
1013 : } else { // if k isn't an integer!
1014 0 : t0 = ( (1-exp(Complex(0, -1.0*C::_2pi*(k-i))))/
1015 0 : (1-exp(Complex(0, -1.0*C::_2pi*(k-i)/Double(ni)))) );
1016 0 : t0 /= ni;
1017 : }
1018 0 : for (size_t j=0; j!=nj; j++) {
1019 0 : if (l_flag) {
1020 0 : t1 = Complex(j==l);
1021 : } else { // if l isn't an integer!
1022 0 : t1 = ( (1-exp(Complex(0, -1.0*C::_2pi*(l-j))))/
1023 0 : (1-exp(Complex(0, -1.0*C::_2pi*(l-j)/Double(nj)))) );
1024 0 : t1 /= nj;
1025 : }
1026 0 : model(i, j) = exp(Complex(0, -1.0*C::_2pi*offset*l_del))*t0*t1;
1027 : }
1028 : }
1029 0 : Complex t2 = sum(data*model);
1030 0 : return t2;
1031 0 : }
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(C::pi*d0)/C::pi)*exp(Complex(0, C::pi*d0));
1056 0 : Complex t1 = Complex(sin(C::pi*d1)/C::pi)*exp(Complex(0, C::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 :
1080 219 : SDBListGridManagerConcat::SDBListGridManagerConcat(SDBList& sdbs) :
1081 219 : SDBListGridManager(sdbs)
1082 : {
1083 219 : std::set<Double> fmaxes;
1084 219 : std::set<Double> fmins;
1085 219 : Float dfn(0.0);
1086 219 : Int totalChans0(0) ;
1087 219 : Int nchan(0);
1088 :
1089 219 : if (sdbs_.nSDB()==0) {
1090 : // The for loop is fine with an empty list, but code below it
1091 : // isn't and there's nothing to lose by bailing early!
1092 0 : return;
1093 : }
1094 :
1095 26754 : for (Int i=0; i != sdbs_.nSDB(); i++) {
1096 26535 : SolveDataBuffer& sdb = sdbs_(i);
1097 26535 : Int spw = sdb.spectralWindow()(0);
1098 26535 : Double t = sdbs_(i).time()(0);
1099 26535 : times_.insert(t);
1100 26535 : if (spwins_.find(spw) == spwins_.end()) {
1101 225 : spwins_.insert(spw);
1102 225 : const Vector<Double>& fs = sdb.freqs();
1103 225 : spwIdToFreqMap_[spw] = &(sdb.freqs());
1104 225 : nchan = sdb.nChannels();
1105 225 : fmaxes.insert(fs(nchan-1));
1106 225 : fmins.insert(fs(0));
1107 : // We assume they're all at the same time.
1108 :
1109 225 : totalChans0 += nchan;
1110 225 : Float df0 = fs(1) - fs(0);
1111 225 : dfn = (fs(nchan-1) - fs(0))/(nchan-1);
1112 : if (DEVDEBUG) {
1113 : cerr << "Spectral window " << spw << " has " << nchan << " channels" << endl;
1114 : cerr << "df0 "<< df0 << "; " << "dfn " << dfn << endl;
1115 : }
1116 : } else {
1117 26310 : continue;
1118 : }
1119 : }
1120 : // nt_ = sdbs_.nSDB()/spwins_.size();
1121 219 : nt_ = times_.size();
1122 219 : tmin_ = *(times_.begin());
1123 219 : tmax_ = *(times_.rbegin());
1124 219 : dt_ = (tmax_ - tmin_)/(nt_ - 1);
1125 219 : nSPWChan_ = nchan;
1126 219 : if (nchan == 1) {
1127 0 : nSPWChan_ = nchan;
1128 0 : fmin_ = *(fmins.begin());
1129 0 : fmax_ = fmin_;
1130 0 : totalChans_ = 1;
1131 0 : df_ = 1;
1132 0 : return;
1133 : }
1134 :
1135 219 : fmin_ = *(fmins.begin());
1136 219 : fmax_ = *(fmaxes.rbegin());
1137 219 : totalChans_ = round((fmax_ - fmin_)/dfn + 1);
1138 219 : df_ = (fmax_ - fmin_)/(totalChans_-1);
1139 : if (DEVDEBUG) {
1140 : cerr << "Global fmin " << fmin_ << " global max " << fmax_ << endl;
1141 : cerr << "nt " << nt_ << " dt " << dt_ << endl;
1142 : cerr << "tmin " << tmin_ << " tmax " << tmax_ << endl;
1143 : cerr << "Global dt " << tmax_ - tmin_ << endl;
1144 : cerr << "Global df " << df_ << endl;
1145 : cerr << "I guess we'll need " << totalChans_ << " freq points in total." << endl;
1146 : cerr << "Compared to " << totalChans0 << " with simple-minded concatenation." << endl;
1147 : cerr << "spwins_.size() " << spwins_.size() << endl;
1148 : cerr << "nSPW() " << nSPW() << endl;
1149 : }
1150 219 : }
1151 :
1152 : // checkAllGridPoints is a diagnostic funtion that should not be called
1153 : // in production releases, but it doesn't do any harm to have it
1154 : // latent.
1155 : void
1156 0 : SDBListGridManagerConcat::checkAllGridpoints() {
1157 0 : map<Int, Vector<Double> const *>::iterator it;
1158 0 : for (it = spwIdToFreqMap_.begin(); it != spwIdToFreqMap_.end(); it++) {
1159 0 : Int spwid = it->first;
1160 0 : Vector<Double> const* fs = it->second;
1161 : Int length;
1162 0 : fs->shape(length);
1163 0 : for (Int i=0; i!=length; i++) {
1164 0 : Double f = (*fs)(i);
1165 0 : Int j = bigFreqGridIndex(f);
1166 : if (DEVDEBUG) {
1167 : cerr << "spwid, i = (" << spwid << ", " << i << ") => " << j << " (" << f << ")" << endl;
1168 : }
1169 : }
1170 : }
1171 : if (DEVDEBUG) {
1172 : cerr << "[1] spwins.size() " << nSPW() << endl;
1173 : }
1174 0 : }
1175 :
1176 : Int
1177 112487 : SDBListGridManagerConcat::swStartIndex(Int spw) {
1178 112487 : Vector<Double> const* fs = spwIdToFreqMap_[spw];
1179 112487 : Double f0 = (*fs)(0);
1180 112487 : return bigFreqGridIndex(f0);
1181 : }
1182 :
1183 219 : DelayRateFFTConcat::DelayRateFFTConcat(SDBList& sdbs, Int refant, Array<Double>& delayWindow,
1184 219 : Array<Double>& rateWindow) :
1185 : DelayRateFFT(sdbs, refant, delayWindow, rateWindow),
1186 219 : gm_(sdbs)
1187 : {
1188 : if (DEVDEBUG) {
1189 : gm_.describe();
1190 : cerr << "gm_.nSPW() " << gm_.nSPW() << endl;
1191 : }
1192 219 : nPadFactor_= max(2, 8 / gm_.nSPW());
1193 219 : nt_ = gm_.nt_;
1194 219 : nPadT_ = nPadFactor_ * nt_;
1195 219 : nChan_ = gm_.nChannels();
1196 219 : nPadChan_ = nPadFactor_*nChan_;
1197 219 : dt_ = gm_.dt_;
1198 219 : f0_ = sdbs.centroidFreq() / 1.e9; // GHz, for delayrate calc
1199 219 : df_ = gm_.df_ / 1.e9;
1200 219 : df_all_ = gm_.fmax_ - gm_.fmin_;
1201 : // This check should be commented out in production:
1202 : // gm_.checkAllGridpoints();
1203 219 : if (nt_ < 2) {
1204 0 : throw(AipsError("Can't do a 2-dimensional FFT on a single timestep! Please consider changing solint to avoid orphan timesteps."));
1205 : }
1206 219 : IPosition ds = delayWindow_.shape();
1207 219 : if (ds.size()!=1 || ds.nelements() != 1) {
1208 0 : throw AipsError("delaywindow must be a list of length 2.");
1209 : }
1210 219 : IPosition rs = rateWindow_.shape();
1211 219 : if (rs.size()!=1 || rs.nelements() != 1) {
1212 0 : throw AipsError("ratewindow must be a list of length 2.");
1213 : }
1214 219 : Int nCorrOrig(sdbs(0).nCorrelations());
1215 219 : nCorr_ = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
1216 575 : for (Int i=0; i<nCorr_; i++) {
1217 356 : activeAntennas_[i].insert(refant_);
1218 : }
1219 :
1220 : // when we get the visCubecorrected it is already
1221 : // reduced to parallel hands, but there isn't a
1222 : // corresponding method for flags.
1223 219 : Int corrStep = (nCorrOrig > 2 ? 3 : 1); // step for p-hands
1224 26754 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
1225 26535 : SolveDataBuffer& s(sdbs(ibuf));
1226 495325 : for (Int irow=0; irow!=s.nRows(); irow++) {
1227 468790 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
1228 468790 : allActiveAntennas_.insert(a1);
1229 468790 : allActiveAntennas_.insert(a2);
1230 : }
1231 : }
1232 219 : nElem_ = 1 + *(allActiveAntennas_.rbegin()) ;
1233 219 : IPosition aggregateDim(2, nCorr_, nElem_);
1234 219 : xcount_.resize(aggregateDim);
1235 219 : sumw_.resize(aggregateDim);
1236 219 : sumww_.resize(aggregateDim);
1237 :
1238 219 : xcount_ = 0;
1239 219 : sumw_ = 0.0;
1240 219 : sumww_ = 0.0;
1241 :
1242 : if (DEVDEBUG) {
1243 : cerr << "Filling FFT grid with " << sdbs.nSDB() << " data buffers." << endl;
1244 : }
1245 219 : IPosition paddedDataSize(4, nCorr_, nElem_, nPadT_, nPadChan_);
1246 219 : Vpad_.resize(paddedDataSize);
1247 219 : Int totalRows = 0;
1248 219 : Int goodRows = 0;
1249 :
1250 :
1251 26754 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
1252 26535 : SolveDataBuffer& s(sdbs(ibuf));
1253 26535 : totalRows += s.nRows();
1254 26535 : if (!s.Ok())
1255 32 : continue;
1256 :
1257 26503 : Int nr = 0;
1258 495101 : for (Int irow=0; irow!=s.nRows(); irow++) {
1259 468598 : if (s.flagRow()(irow))
1260 356111 : continue;
1261 : Int iant;
1262 466945 : Int a1(s.antenna1()(irow)), a2(s.antenna2()(irow));
1263 466945 : if (a1 == a2) {
1264 52017 : continue;
1265 : }
1266 414928 : else if (a1 == refant_) {
1267 112487 : iant = a2;
1268 : }
1269 302441 : else if (a2 == refant_) {
1270 0 : iant = a1;
1271 : }
1272 : else {
1273 302441 : continue;
1274 : }
1275 : // OK, we're not skipping this one so we have to do something.
1276 :
1277 : // v has shape (nelems, ?, nrows, nchannels)
1278 112487 : Cube<Complex> v = s.visCubeCorrected();
1279 112487 : Cube<Float> w = s.weightSpectrum();
1280 112487 : Cube<Bool> fl = s.flagCube();
1281 112487 : Int spw = s.spectralWindow()(0);
1282 112487 : Int f_index = gm_.swStartIndex(spw); // ditto!
1283 112487 : Int t_index = gm_.getTimeIndex(s.time()(0));
1284 112487 : Int spwchans = gm_.nSPWChan_;
1285 112487 : IPosition start(4, 0, iant, t_index, f_index);
1286 112487 : IPosition stop(4, nCorr_, 1, 1, spwchans);
1287 112487 : IPosition stride(4, 1, 1, 1, 1);
1288 112487 : Slicer sl1(start, stop, stride, Slicer::endIsLength);
1289 224974 : Slicer sl2(IPosition(3, 0, 0, irow),
1290 224974 : IPosition(3, nCorr_, spwchans, 1),
1291 224974 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
1292 :
1293 224974 : Slicer flagSlice(IPosition(3, 0, 0, irow),
1294 224974 : IPosition(3, nCorr_, spwchans, 1),
1295 224974 : IPosition(3, corrStep, 1, 1), Slicer::endIsLength);
1296 112487 : nr++;
1297 : if (DEVDEBUG && false) {
1298 : cerr << "nr " << nr
1299 : << " irow " << endl
1300 : << "Vpad shape " << Vpad_.shape() << endl
1301 : << "v shape " << v.shape() << endl
1302 : << "sl2 " << sl2 << endl
1303 : << "sl1 " << sl1 << endl
1304 : << "flagSlice " << flagSlice << endl;
1305 : }
1306 112487 : Array<Complex> rhs = v(sl2).nonDegenerate(1);
1307 112487 : Array<Float> weights = w(sl2).nonDegenerate(1);
1308 :
1309 112487 : unitize(rhs);
1310 112487 : Vpad_(sl1).nonDegenerate(1) = rhs * weights;
1311 :
1312 112487 : Array<Bool> flagged(fl(flagSlice).nonDegenerate(1));
1313 : // Zero flagged entries.
1314 112487 : Vpad_(sl1).nonDegenerate(1)(flagged) = Complex(0.0);
1315 :
1316 112487 : if (!allTrue(flagged)) {
1317 290971 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
1318 178484 : IPosition p(2, icorr, iant);
1319 178484 : Bool actually = false;
1320 178484 : activeAntennas_[icorr].insert(iant);
1321 3014932 : for (Int ichan=0; ichan != (Int) spwchans; ichan++) {
1322 2836448 : IPosition pchan(2, icorr, ichan);
1323 2836448 : if (!flagged(pchan)) {
1324 2795500 : Float wv = weights(pchan);
1325 2795500 : if (wv < 0) {
1326 0 : cerr << "spwchans " << spwchans << endl;
1327 0 : cerr << "Negative weight << (" << wv << ") on row "
1328 0 : << irow << " baseline (" << a1 << ", " << a2 << ") "
1329 0 : << " channel " << ichan << endl;
1330 0 : cerr << "pchan " << pchan << endl;
1331 0 : cerr << "Weights " << weights << endl;
1332 : }
1333 2795500 : xcount_(p)++;
1334 2795500 : sumw_(p) += wv;
1335 2795500 : sumww_(p) += wv*wv;
1336 2795500 : actually = true;
1337 : }
1338 2836448 : }
1339 178484 : if (actually) {
1340 178484 : activeAntennas_[icorr].insert(iant);
1341 178484 : goodRows++;
1342 : }
1343 178484 : }
1344 : }
1345 : if (DEVDEBUG && 0) {
1346 : cerr << "flagged " << flagged << endl;
1347 : cerr << "flagSlice " << flagSlice << endl
1348 : << "fl.shape() " << fl.shape() << endl
1349 : << "Vpad_.shape() " << Vpad_.shape() << endl
1350 : << "flagged.shape() " << flagged.shape() << endl
1351 : << "sl1 " << sl1 << endl;
1352 : }
1353 112487 : }
1354 : }
1355 : if (DEVDEBUG) {
1356 : cerr << "In DelayRateFFTConcat::DelayRateFFTConcat " << endl;
1357 : printActive();
1358 : cerr << "sumw_ " << sumw_ << endl;
1359 : cerr << "Constructed a DelayRateFFTConcat object." << endl;
1360 : cerr << "totalRows " << totalRows << endl;
1361 : cerr << "goodRows " << goodRows << endl;
1362 : }
1363 :
1364 219 : }
1365 :
1366 0 : DelayRateFFTConcat::DelayRateFFTConcat(Array<Complex>& data, Int nPadFactor, Float f0, Float df, Float dt, SDBList& s,
1367 0 : Array<Double>& delayWindow, Array<Double>& rateWindow) :
1368 : DelayRateFFT(s, 0, delayWindow, rateWindow),
1369 0 : gm_(s),
1370 0 : nPadFactor_(nPadFactor),
1371 0 : Vpad_() {
1372 0 : dt_ = dt;
1373 0 : f0_ = f0;
1374 0 : df_ = df;
1375 0 : IPosition shape = data.shape();
1376 0 : nCorr_ = shape(0);
1377 0 : nElem_ = shape(1);
1378 0 : nt_ = shape(2);
1379 0 : nChan_ = shape(3);
1380 0 : nPadT_ = nPadFactor_*nt_;
1381 0 : nPadChan_ = nPadFactor_*nChan_;
1382 0 : IPosition paddedDataSize(4, nCorr_, nElem_, nPadT_, nPadChan_);
1383 0 : Vpad_.resize(paddedDataSize);
1384 :
1385 0 : IPosition start(4, 0, 0, 0, 0);
1386 0 : IPosition stop(4, nCorr_, nElem_, nt_, nChan_);
1387 0 : IPosition stride(4, 1, 1, 1, 1);
1388 0 : Slicer sl1(start, stop, stride, Slicer::endIsLength);
1389 0 : Vpad_(sl1) = data;
1390 :
1391 0 : unitize(Vpad_);
1392 :
1393 0 : }
1394 :
1395 : void
1396 219 : DelayRateFFTConcat::FFT() {
1397 : if (DEVDEBUG) {
1398 : cerr << "DelayRateFFTConcat::FFT()" << endl;
1399 : }
1400 : // Axes are 0: correlation (i.e., hand of polarization), 1: antenna, 2: time, 3: channel
1401 219 : Vector<Bool> ax(4, false);
1402 219 : ax(2) = true;
1403 219 : ax(3) = true;
1404 : // Also copied from DelayFFT in KJones.
1405 : // we make a copy to FFT in place.
1406 : if (DEVDEBUG) {
1407 : cerr << "Vpad_.shape() " << Vpad_.shape() << endl;
1408 : }
1409 219 : ArrayLattice<Complex> c(Vpad_);
1410 219 : LatticeFFT::cfft0(c, ax, true);
1411 : if (DEVDEBUG) {
1412 : cerr << "FFT transformed" << endl;
1413 : }
1414 219 : }
1415 :
1416 : std::pair<Bool, Float>
1417 4088 : DelayRateFFTConcat::xinterp(Float alo, Float amax, Float ahi) {
1418 4088 : if (amax > 0.0 && alo == amax && amax == ahi)
1419 0 : return std::make_pair(true, 0.5);
1420 :
1421 4088 : Float denom(alo-2.0*amax+ahi);
1422 4088 : Bool cond = amax>0.0 && abs(denom)>0.0 ;
1423 4088 : Float fpk = cond ? 0.5-(ahi-amax)/denom : 0.0;
1424 4088 : return std::make_pair(cond, fpk);
1425 : }
1426 :
1427 : void
1428 219 : DelayRateFFTConcat::searchPeak() {
1429 : if (DEVDEBUG) {
1430 : cerr << "DelayRateFFTConcat::searchPeak()" << endl;
1431 : }
1432 :
1433 : // Recall param_ -> [phase, delay, rate] for each correlation
1434 219 : param_.resize(3*nCorr_, nElem_); // This might be better done elsewhere.
1435 219 : param_.set(0.0);
1436 219 : flag_.resize(3*nCorr_, nElem_);
1437 219 : flag_.set(true); // all flagged initially
1438 : if (DEVDEBUG) {
1439 : cerr << "nt_ " << nt_ << " nPadChan_ " << nPadChan_ << endl;
1440 : cerr << "Vpad_.shape() " << Vpad_.shape() << endl;
1441 : cerr << "delayWindow_ " << delayWindow_ << endl;
1442 :
1443 : }
1444 :
1445 575 : for (Int icorr=0; icorr<nCorr_; ++icorr) {
1446 356 : flag_(icorr*3 + 0, refant()) = false;
1447 356 : flag_(icorr*3 + 1, refant()) = false;
1448 356 : flag_(icorr*3 + 2, refant()) = false;
1449 2756 : for (Int ielem=0; ielem<nElem_; ++ielem) {
1450 2400 : if (ielem==refant()) {
1451 356 : continue;
1452 : }
1453 : // NB: Time, Channel
1454 : // And once again we fail at slicing
1455 2044 : IPosition start(4, icorr, ielem, 0, 0);
1456 2044 : IPosition stop(4, 1, 1, nPadT_, nPadChan_);
1457 2044 : IPosition step(4, 1, 1, 1, 1);
1458 2044 : Slicer sl(start, stop, step, Slicer::endIsLength);
1459 4088 : Matrix<Complex> aS = Vpad_(sl).nonDegenerate();
1460 2044 : Int sgn = (ielem < refant()) ? 1 : -1;
1461 :
1462 : // Below is the gory details for turning delay window into index range
1463 2044 : Double bw = Float(nPadChan_)*df_;
1464 2044 : Double d0 = sgn*delayWindow_(IPosition(1, 0));
1465 2044 : Double d1 = sgn*delayWindow_(IPosition(1, 1));
1466 2044 : if (d0 > d1) std::swap(d0, d1);
1467 2044 : d0 = max(d0, -0.5/df_);
1468 2044 : d1 = min(d1, (0.5-1/nPadChan_)/df_);
1469 :
1470 : // It's simpler to keep the ranges as signed integers and
1471 : // handle the wrapping of the FFT in the loop over
1472 : // indices. Recall that the FFT result returned has indices
1473 : // that run from 0 to nPadChan_/2 -1 and then from
1474 : // -nPadChan/2 to -1, so far as our delay is concerned.
1475 2044 : Int i0 = bw*d0;
1476 2044 : Int i1 = bw*d1;
1477 2044 : if (i1==i0) i1++;
1478 : // Now for the gory details for turning rate window into index range
1479 2044 : Double width = nPadT_*dt_*1e9*f0_;
1480 2044 : Double r0 = sgn*rateWindow_(IPosition(1,0));
1481 2044 : Double r1 = sgn*rateWindow_(IPosition(1,1));
1482 2044 : if (r0 > r1) std::swap(r0, r1);
1483 2044 : r0 = max(r0, -0.5/(dt_*1e9*f0_));
1484 2044 : r1 = min(r1, (0.5 - 1/nPadT_)/(dt_*1e9*f0_));
1485 :
1486 2044 : Int j0 = width*r0;
1487 2044 : Int j1 = width*r1;
1488 2044 : if (j1==j0) j1++;
1489 : if (DEVDEBUG) {
1490 : cerr << "Checking the windows for delay and rate search." << endl;
1491 : cerr << "bw " << bw << endl;
1492 : cerr << "d0 " << d0 << " d1 " << d1 << endl;
1493 : cerr << "i0 " << i0 << " i1 " << i1 << endl;
1494 : cerr << "r0 " << r0 << " r1 " << r1 << endl;
1495 : cerr << "j0 " << j0 << " j1 " << j1 << endl;
1496 : }
1497 2044 : Matrix<Float> amp(amplitude(aS));
1498 2044 : Int ipkch(0);
1499 2044 : Int ipkt(0);
1500 2044 : Float amax(-1.0);
1501 : // Unlike KJones we have to iterate in time too
1502 1670562 : for (Int itime0=j0; itime0 != j1; itime0++) {
1503 1668518 : Int itime = (itime0 < 0) ? itime0 + nPadT_ : itime0;
1504 233652006 : for (Int ich0=i0; ich0 != i1; ich0++) {
1505 231983488 : Int ich = (ich0 < 0) ? ich0 + nPadChan_ : ich0;
1506 : // cerr << "Gridpoint " << itime << ", " << ich << "->" << amp(itime, ich) << endl;
1507 231983488 : if (amp(itime, ich) > amax) {
1508 197509 : ipkch = ich;
1509 197509 : ipkt = itime;
1510 197509 : amax=amp(itime, ich);
1511 : }
1512 : }
1513 : }
1514 : // Finished grovelling. Now we have the location of the
1515 : // maximum amplitude.
1516 2044 : Float alo_ch = amp(ipkt, (ipkch > 0) ? ipkch-1 : nPadChan_-1);
1517 2044 : Float ahi_ch = amp(ipkt, ipkch<(nPadChan_-1) ? ipkch+1 : 0);
1518 2044 : std::pair<Bool, Float> maybeFpkch = xinterp(alo_ch, amax, ahi_ch);
1519 : // We handle wrapping while looking for neighbours
1520 2044 : Float alo_t = amp(ipkt > 0 ? ipkt-1 : nPadT_ -1, ipkch);
1521 2044 : Float ahi_t = amp(ipkt < (nPadT_ -1) ? ipkt+1 : 0, ipkch);
1522 : if (DEVDEBUG) {
1523 : cerr << "For element " << ielem << endl;
1524 : cerr << "In channel dimension ipkch " << ipkch << " alo " << alo_ch
1525 : << " amax " << amax << " ahi " << ahi_ch << endl;
1526 : cerr << "In time dimension ipkt " << ipkt << " alo " << alo_t
1527 : << " amax " << amax << " ahi " << ahi_t << endl;
1528 : }
1529 2044 : std::pair<Bool, Float> maybeFpkt = xinterp(alo_t, amax, ahi_t);
1530 :
1531 2044 : if (maybeFpkch.first and maybeFpkt.first) {
1532 : // Phase
1533 1800 : Complex c = aS(ipkt, ipkch);
1534 1800 : Float phase = arg(c);
1535 1800 : param_(icorr*3 + 0, ielem) = sgn*phase;
1536 1800 : Float delay = (ipkch)/Float(nPadChan_);
1537 1800 : if (delay > 0.5) delay -= 1.0; // fold
1538 1800 : delay /= df_; // nsec
1539 1800 : param_(icorr*3 + 1, ielem) = sgn*delay; //
1540 1800 : Double rate = (ipkt)/Float(nPadT_);
1541 1800 : if (rate > 0.5) rate -= 1.0;
1542 1800 : Double rate0 = rate/dt_;
1543 1800 : Double rate1 = rate0/(1e9 * f0_);
1544 :
1545 1800 : param_(icorr*3 + 2, ielem) = Float(sgn*rate1);
1546 : if (DEVDEBUG) {
1547 : cerr << "maybeFpkch.second=" << maybeFpkch.second
1548 : << ", df_= " << df_
1549 : << " fpkch=" << (ipkch + maybeFpkch.second) << endl;
1550 : cerr << " maybeFpkt.second=" << maybeFpkt.second
1551 : << " rate0=" << rate
1552 : << " 1e9 * f0_=" << 1e9 * f0_
1553 : << ", dt_=" << dt_
1554 : << " fpkt=" << (ipkt + maybeFpkt.second) << endl;
1555 :
1556 : }
1557 : if (DEVDEBUG) {
1558 : cerr << "Found peak for element " << ielem << " correlation " << icorr
1559 : << " ipkt=" << ipkt << "/" << nPadT_ << ", ipkch=" << ipkch << "/" << nPadChan_
1560 : << " peak=" << amax
1561 : << "; delay " << delay << ", rate " << rate
1562 : << ", phase " << arg(c) << " sign= " << sgn << endl;
1563 : }
1564 : // Set 3 flags.
1565 1800 : flag_(icorr*3 + 0, ielem)=false;
1566 1800 : flag_(icorr*3 + 1, ielem)=false;
1567 1800 : flag_(icorr*3 + 2, ielem)=false;
1568 1800 : }
1569 : else {
1570 : if (DEVDEBUG) {
1571 : cerr << "No peak in 2D FFT for element " << ielem << " correlation " << icorr << endl;
1572 : }
1573 : }
1574 2044 : }
1575 : }
1576 219 : }
1577 :
1578 :
1579 : Float
1580 1800 : DelayRateFFTConcat::snr(Int icorr, Int ielem, Float delay, Float rate) {
1581 : if (DEVDEBUG) {
1582 : cerr << "DelayRateFFTConcat::snr"<< endl;
1583 : }
1584 : // We calculate a signal-to-noise ration for the 2D FFT fringefit
1585 : // using a formula transcribed from AIPS FRING.
1586 : //
1587 : // Have to convert delay and rate back into indices on the padded 2D grid.
1588 1800 : Int sgn = (ielem < refant()) ? 1 : -1;
1589 1800 : delay *= sgn*df_;
1590 1800 : if (delay < 0.0) delay += 1;
1591 1800 : Int ichan = Int(delay*nPadChan_ + 0.5);
1592 1800 : if (ichan == nPadChan_) ichan = 0;
1593 :
1594 1800 : rate *= sgn*1e9 * f0_;
1595 1800 : rate *= dt_;
1596 1800 : if (rate < 0.0) rate += 1;
1597 1800 : Int itime = Int(rate*nPadT_ + 0.5);
1598 1800 : if (itime == nPadT_) itime = 0;
1599 : // What about flags? If the datapoint closest to the computed
1600 : // delay and rate values is flagged we probably shouldn't use
1601 : // it, but what *should* we use?
1602 1800 : IPosition ipos(4, icorr, ielem, itime, ichan);
1603 1800 : IPosition p(2, icorr, ielem);
1604 1800 : Complex v = Vpad_(ipos);
1605 1800 : Float peak = abs(v);
1606 1800 : if (peak > 0.999*sumw_(p)) peak=0.999*sumw_(p);
1607 : // xcount is number of data points for baseline to ielem
1608 : // sumw is sum of weights,
1609 : // sumww is sum of squares of weights
1610 :
1611 : Float cwt;
1612 1800 : if (fabs(sumw_(p))<FLT_EPSILON) {
1613 0 : cwt = 0;
1614 : if (DEVDEBUG) {
1615 : cerr << "Correlation " << icorr << " antenna " << ielem << ": sum of weights is zero." << endl;
1616 : }
1617 : } else {
1618 1800 : Float x = C::pi/2*peak/sumw_(p);
1619 : // The magic numbers in the following formula are from AIPS FRING
1620 1800 : cwt = (pow(tan(x), 1.163) * sqrt(sumw_(p)/sqrt(sumww_(p)/xcount_(p))));
1621 : if (DEVDEBUG) {
1622 : cerr << "Correlation " << icorr << " antenna " << ielem << " ipos " << ipos
1623 : << " peak=" << peak << "; xang=" << x << "; xcount=" << xcount_(p) << "; sumw=" << sumw_(p) << "; sumww=" << sumww_(p)
1624 : << " snr " << cwt << endl;
1625 : }
1626 : }
1627 1800 : return cwt;
1628 1800 : }
1629 :
1630 :
1631 : } // end namespace
|
