Source
1423
1423
}
1424
1424
}
1425
1425
gsl_matrix_free(lu);
1426
1426
gsl_matrix_free(inv_hess);
1427
1427
// SNR[i], according to aips, is 1/sqrt(sigma2*hess(i1,i1)*0.5);
1428
1428
// Note that in AIPS i1 ranges over 1..NANT
1429
1429
// We use 1 as a success code.
1430
1430
return 1;
1431
1431
}
1432
1432
1433
+
1434
+
Int
1435
+
findRefAntWithData(SDBList& sdbs, Vector<Int>& refAntList, Int prtlev) {
1436
+
std::set<Int> activeAntennas;
1437
+
for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
1438
+
SolveDataBuffer& s(sdbs(ibuf));
1439
+
if (!s.Ok())
1440
+
continue;
1441
+
Cube<Bool> fl = s.flagCube();
1442
+
for (Int irow=0; irow!=s.nRows(); irow++) {
1443
+
if (s.flagRow()(irow))
1444
+
continue;
1445
+
Int a1(s.antenna1()(irow));
1446
+
Int a2(s.antenna2()(irow));
1447
+
// Not using irow
1448
+
Matrix<Bool> flr = fl.xyPlane(irow);
1449
+
if (!allTrue(flr)) {
1450
+
activeAntennas.insert(a1);
1451
+
activeAntennas.insert(a2);
1452
+
}
1453
+
}
1454
+
}
1455
+
if (prtlev > 2) {
1456
+
cout << "[FringeJones.cc::findRefAntWithData] activeAntennas: ";
1457
+
std::copy(
1458
+
activeAntennas.begin(),
1459
+
activeAntennas.end(),
1460
+
std::ostream_iterator<Int>(std::cout, " ")
1461
+
);
1462
+
cout << endl;
1463
+
}
1464
+
Int refAnt = -1;
1465
+
for (Vector<Int>::ConstIteratorSTL a = refAntList.begin(); a != refAntList.end(); a++) {
1466
+
if (activeAntennas.find(*a) != activeAntennas.end()) {
1467
+
if (prtlev > 2)
1468
+
cout << "[FringeJones.cc::findRefAntWithData] We are choosing refant " << *a << endl;
1469
+
refAnt = *a;
1470
+
break;
1471
+
} else {
1472
+
if (prtlev > 2)
1473
+
cout << "[FringeJones.cc::findRefAntWithData] No data for refant " << *a << endl;
1474
+
}
1475
+
}
1476
+
return refAntList(0);
1477
+
}
1478
+
1479
+
1433
1480
// Stolen from SolveDataBuffer
1434
1481
void
1435
1482
aggregateTimeCentroid(SDBList& sdbs, Int refAnt, std::map<Int, Double>& aggregateTime) {
1436
1483
// Weighted average of SDBs' timeCentroids
1437
1484
std::map<Int, Double> aggregateWeight;
1438
1485
for (Int i=0; i < sdbs.nSDB(); ++i) {
1439
1486
SolveDataBuffer& s = sdbs(i);
1440
1487
Vector<Double> wtvD;
1441
1488
// Sum over correlations and channels to get a vector of weights for each row
1442
1489
Vector<Float> wtv(partialSums(s.weightSpectrum(), IPosition(2, 0, 1)));
1886
1933
for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1887
1934
if (spwMap()(ispw)>-1)
1888
1935
KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1889
1936
}
1890
1937
1891
1938
}
1892
1939
1893
1940
void FringeJones::setSolve(const Record& solve) {
1894
1941
1895
1942
// Call parent to do conventional things
1943
+
if (prtlev() > 2) {
1944
+
cout << "Before GJones::setSolve" << endl
1945
+
<< "FringeJones::setSolve()" <<endl
1946
+
<< "FringeJones::refant() = "<< refant() <<endl
1947
+
<< "FringeJones::refantlist() = "<< refantlist() <<endl;
1948
+
}
1896
1949
GJones::setSolve(solve);
1897
1950
1898
1951
// if (!ct_)
1899
1952
// throw(AipsError("No calibration table specified"));
1900
1953
// cerr << "setSolve here, ct_: "<< ct_ << endl;
1901
1954
1902
1955
// Trap unspecified refant:
1903
1956
if (refant()<0)
1904
1957
throw(AipsError("Please specify a good reference antenna (refant) explicitly."));
1905
1958
if (solve.isDefined("zerorates")) {
2037
2090
Double t0 = sdbs(0).time()(0);
2038
2091
Double dt0 = refTime() - t0;
2039
2092
//Double df0 = ref_freq - sdbs.freqs()(0);
2040
2093
//Double df0 = 0;
2041
2094
Double df0 = centroidFreq - sdbs(0).freqs()(0); // global center to global edge
2042
2095
2043
2096
logSink() << "Solving for fringes for spw=" << currSpw() << " at t="
2044
2097
<< MVTime(refTime()/C::day).string(MVTime::YMD,7) << LogIO::POST;
2045
2098
2046
2099
std::map<Int, Double> aggregateTime;
2100
+
// Set the refant to the first choice that has data!
2101
+
refant() = findRefAntWithData(sdbs, refantlist(), prtlev());
2102
+
if (refant()<0)
2103
+
throw(AipsError("No valid reference antenna supplied."));
2104
+
2047
2105
aggregateTimeCentroid(sdbs, refant(), aggregateTime);
2048
2106
2049
2107
if (DEVDEBUG) {
2050
2108
std::cerr << "Weighted time centroids" << endl;
2051
2109
for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it)
2052
2110
std::cerr << it->first << " => " << it->second - t0 << std::endl;
2053
2111
}
2054
2112
2055
2113
DelayRateFFT drf(sdbs, refant(), delayWindow(), rateWindow());
2056
2114
drf.FFT();