Source
22
22
//# National Radio Astronomy Observatory
23
23
//# 520 Edgemont Road
24
24
//# Charlottesville, VA 22903-2475 USA
25
25
//#
26
26
27
27
#include <synthesis/MeasurementComponents/FringeJones.h>
28
28
29
29
#include <msvis/MSVis/VisBuffer.h>
30
30
#include <msvis/MSVis/VisBuffAccumulator.h>
31
31
#include <ms/MeasurementSets/MSColumns.h>
32
+
#include <synthesis/CalTables/CTIter.h>
32
33
#include <synthesis/MeasurementEquations/VisEquation.h> // *
33
34
#include <synthesis/MeasurementComponents/SolveDataBuffer.h>
35
+
#include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
34
36
#include <lattices/Lattices/ArrayLattice.h>
35
37
#include <lattices/LatticeMath/LatticeFFT.h>
36
38
#include <scimath/Mathematics/FFTServer.h>
37
39
38
40
#include <casa/Arrays/ArrayMath.h>
39
41
#include <casa/Arrays/MatrixMath.h>
40
42
#include <casa/Arrays/ArrayLogical.h>
41
43
#include <casa/BasicSL/String.h>
42
44
#include <casa/Utilities/Assert.h>
43
45
#include <casa/Exceptions/Error.h>
2178
2180
}
2179
2181
}
2180
2182
}
2181
2183
2182
2184
void
2183
2185
FringeJones::solveOneVB(const VisBuffer&) {
2184
2186
throw(AipsError("VisBuffer interface not supported!"));
2185
2187
}
2186
2188
2187
2189
2190
+
void FringeJones::globalPostSolveTinker() {
2191
+
2192
+
// Re-reference the phase, if requested
2193
+
if (refantlist()(0)>-1) applyRefAnt();
2194
+
}
2195
+
2196
+
void FringeJones::applyRefAnt() {
2197
+
2198
+
// TBD:
2199
+
// 1. Synchronize refant changes on par axis
2200
+
// 2. Implement minimum mean deviation algorithm
2201
+
2202
+
if (refantlist()(0)<0)
2203
+
throw(AipsError("No refant specified."));
2204
+
2205
+
Int nUserRefant=refantlist().nelements();
2206
+
2207
+
// Get the preferred refant names from the MS
2208
+
String refantName(msmc().antennaName(refantlist()(0)));
2209
+
if (nUserRefant>1) {
2210
+
refantName+=" (";
2211
+
for (Int i=1;i<nUserRefant;++i) {
2212
+
refantName+=msmc().antennaName(refantlist()(i));
2213
+
if (i<nUserRefant-1) refantName+=",";
2214
+
}
2215
+
refantName+=")";
2216
+
}
2217
+
2218
+
logSink() << "Applying refant: " << refantName
2219
+
<< " refantmode = " << refantmode();
2220
+
if (refantmode()=="flex")
2221
+
logSink() << " (hold alternate refants' phase constant) when refant flagged";
2222
+
if (refantmode()=="strict")
2223
+
logSink() << " (flag all antennas when refant flagged)";
2224
+
logSink() << LogIO::POST;
2225
+
2226
+
// Generate a prioritized refant choice list
2227
+
// The first entry in this list is the user's primary refant,
2228
+
// the second entry is the refant used on the previous interval,
2229
+
// and the rest is a prioritized list of alternate refants,
2230
+
// starting with the user's secondary (if provided) refants,
2231
+
// followed by the rest of the array, in distance order. This
2232
+
// makes the priorities correct all the time, and prevents
2233
+
// a semi-stochastic alternation (by preferring the last-used
2234
+
// alternate, even if nominally higher-priority refants become
2235
+
// available)
2236
+
2237
+
2238
+
// Extract antenna positions
2239
+
Matrix<Double> xyz;
2240
+
if (msName()!="<noms>") {
2241
+
MeasurementSet ms(msName());
2242
+
ROMSAntennaColumns msant(ms.antenna());
2243
+
msant.position().getColumn(xyz);
2244
+
}
2245
+
else {
2246
+
// TBD RO*
2247
+
CTColumns ctcol(*ct_);
2248
+
CTAntennaColumns& antcol(ctcol.antenna());
2249
+
antcol.position().getColumn(xyz);
2250
+
}
2251
+
2252
+
// Calculate (squared) antenna distances, relative
2253
+
// to last preferred antenna
2254
+
Vector<Double> dist2(xyz.ncolumn(),0.0);
2255
+
for (Int i=0;i<3;++i) {
2256
+
Vector<Double> row=xyz.row(i);
2257
+
row-=row(refantlist()(nUserRefant-1));
2258
+
dist2+=square(row);
2259
+
}
2260
+
// Move preferred antennas to a large distance
2261
+
for (Int i=0;i<nUserRefant;++i)
2262
+
dist2(refantlist()(i))=DBL_MAX;
2263
+
2264
+
// Generated sorted index
2265
+
Vector<uInt> ord;
2266
+
genSort(ord,dist2);
2267
+
2268
+
// Assemble the whole choices list
2269
+
Int nchoices=nUserRefant+1+ord.nelements();
2270
+
Vector<Int> refantchoices(nchoices,0);
2271
+
Vector<Int> r(refantchoices(IPosition(1,nUserRefant+1),IPosition(1,refantchoices.nelements()-1)));
2272
+
convertArray(r,ord);
2273
+
2274
+
// set first two to primary preferred refant
2275
+
refantchoices(0)=refantchoices(1)=refantlist()(0);
2276
+
2277
+
// set user's secondary refants (if any)
2278
+
if (nUserRefant>1)
2279
+
refantchoices(IPosition(1,2),IPosition(1,nUserRefant))=
2280
+
refantlist()(IPosition(1,1),IPosition(1,nUserRefant-1));
2281
+
2282
+
//cout << "refantchoices = " << refantchoices << endl;
2283
+
2284
+
2285
+
if (refantmode()=="strict") {
2286
+
nchoices=1;
2287
+
refantchoices.resize(1,True);
2288
+
}
2289
+
2290
+
Vector<Int> nPol(nSpw(),nPar()); // TBD:or 1, if data was single pol
2291
+
2292
+
if (nPar()==6) {
2293
+
// Verify that 2nd poln has unflagged solutions, PER SPW
2294
+
ROCTMainColumns ctmc(*ct_);
2295
+
2296
+
Block<String> cols(1);
2297
+
cols[0]="SPECTRAL_WINDOW_ID";
2298
+
CTIter ctiter(*ct_,cols);
2299
+
Cube<Bool> fl;
2300
+
2301
+
while (!ctiter.pastEnd()) {
2302
+
2303
+
Int ispw=ctiter.thisSpw();
2304
+
fl.assign(ctiter.flag());
2305
+
2306
+
IPosition blc(3,0,0,0), trc(fl.shape());
2307
+
trc-=1; trc(0)=blc(0)=1;
2308
+
2309
+
// cout << "ispw = " << ispw << " nfalse(fl(1,:,:)) = " << nfalse(fl(blc,trc)) << endl;
2310
+
2311
+
// If there are no unflagged solutions in 2nd pol,
2312
+
// avoid it in refant calculations
2313
+
if (nfalse(fl(blc,trc))==0)
2314
+
nPol(ispw)=1;
2315
+
2316
+
ctiter.next();
2317
+
}
2318
+
}
2319
+
// cout << "nPol = " << nPol << endl;
2320
+
2321
+
Bool usedaltrefant(false);
2322
+
Int currrefant(refantchoices(0)), lastrefant(-1);
2323
+
2324
+
Block<String> cols(2);
2325
+
cols[0]="SPECTRAL_WINDOW_ID";
2326
+
cols[1]="TIME";
2327
+
CTIter ctiter(*ct_,cols);
2328
+
2329
+
// Arrays to hold per-timestamp solutions
2330
+
Cube<Float> solA, solB;
2331
+
Cube<Bool> flA, flB;
2332
+
Vector<Int> ant1A, ant1B, ant2B;
2333
+
Matrix<Complex> refPhsr; // the reference phasor [npol,nchan]
2334
+
Int lastspw(-1);
2335
+
Bool first(true);
2336
+
while (!ctiter.pastEnd()) {
2337
+
Int ispw=ctiter.thisSpw();
2338
+
if (ispw!=lastspw) first=true; // spw changed, start over
2339
+
2340
+
// Read in the current sol, fl, ant1:
2341
+
solB.assign(ctiter.fparam());
2342
+
flB.assign(ctiter.flag());
2343
+
ant1B.assign(ctiter.antenna1());
2344
+
ant2B.assign(ctiter.antenna2());
2345
+
2346
+
// First time thru, 'previous' solution same as 'current'
2347
+
if (first) {
2348
+
solA.reference(solB);
2349
+
flA.reference(flB);
2350
+
ant1A.reference(ant1B);
2351
+
}
2352
+
IPosition shB(solB.shape());
2353
+
IPosition shA(solA.shape());
2354
+
2355
+
// Find a good refant at this time
2356
+
// A good refant is one that is unflagged in all polarizations
2357
+
// in the current(B) and previous(A) intervals (so they can be connected)
2358
+
Int irefA(0),irefB(0); // index on 3rd axis of solution arrays
2359
+
Int ichoice(0); // index in refantchoicelist
2360
+
Bool found(false);
2361
+
IPosition blcA(3,0,0,0),trcA(shA),blcB(3,0,0,0),trcB(shB);
2362
+
trcA-=1; trcA(0)=trcA(2)=0;
2363
+
trcB-=1; trcB(0)=trcB(2)=0;
2364
+
ichoice=0;
2365
+
while (!found && ichoice<nchoices) {
2366
+
// Find index of current refant choice
2367
+
irefA=irefB=0;
2368
+
while (ant1A(irefA)!=refantchoices(ichoice) && irefA<shA(2)) ++irefA;
2369
+
while (ant1B(irefB)!=refantchoices(ichoice) && irefB<shB(2)) ++irefB;
2370
+
2371
+
if (irefA<shA(2) && irefB<shB(2)) {
2372
+
2373
+
// cout << " Trial irefA,irefB: " << irefA << "," << irefB
2374
+
// << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2375
+
2376
+
blcA(2)=trcA(2)=irefA;
2377
+
blcB(2)=trcB(2)=irefB;
2378
+
found=true; // maybe
2379
+
for (Int ipol=0;ipol<nPol(ispw);++ipol) {
2380
+
blcA(0)=trcA(0)=blcB(0)=trcB(0)=ipol;
2381
+
found &= (nfalse(flA(blcA,trcA))>0); // previous interval
2382
+
found &= (nfalse(flB(blcB,trcB))>0); // current interval
2383
+
}
2384
+
}
2385
+
else
2386
+
// irefA or irefB out-of-range
2387
+
found=false; // Just to be sure
2388
+
2389
+
if (!found) ++ichoice; // try next choice next round
2390
+
2391
+
}
2392
+
2393
+
if (found) {
2394
+
// at this point, irefA/irefB point to a good refant
2395
+
2396
+
// Keep track
2397
+
usedaltrefant|=(ichoice>0);
2398
+
currrefant=refantchoices(ichoice);
2399
+
refantchoices(1)=currrefant; // 2nd priorty next time
2400
+
2401
+
// cout << " currrefant = " << currrefant << " (" << ichoice << ")" << endl;
2402
+
2403
+
// cout << " Final irefA,irefB: " << irefA << "," << irefB
2404
+
// << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2405
+
2406
+
2407
+
// Only report if using an alternate refant
2408
+
if (currrefant!=lastrefant && ichoice>0) {
2409
+
logSink()
2410
+
<< "At "
2411
+
<< MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2412
+
<< " ("
2413
+
<< "Spw=" << ctiter.thisSpw()
2414
+
<< ", Fld=" << ctiter.thisField()
2415
+
<< ")"
2416
+
<< ", using refant " << msmc().antennaName(currrefant)
2417
+
<< " (id=" << currrefant
2418
+
<< ")" << " (alternate)"
2419
+
<< LogIO::POST;
2420
+
}
2421
+
2422
+
// Form reference phasor [nPar,nChan]
2423
+
Matrix<Float> rA,rB;
2424
+
Matrix<Bool> rflA,rflB;
2425
+
rB.assign(solB.xyPlane(irefB));
2426
+
rflB.assign(flB.xyPlane(irefB));
2427
+
2428
+
if (!first) {
2429
+
// Get and condition previous phasor for the current refant
2430
+
rA.assign(solA.xyPlane(irefA));
2431
+
rflA.assign(flA.xyPlane(irefA));
2432
+
rB-=rA;
2433
+
2434
+
// Accumulate flags
2435
+
rflB&=rflA;
2436
+
}
2437
+
2438
+
// cout << " rB = " << rB << endl;
2439
+
// cout << boolalpha << " rflB = " << rflB << endl;
2440
+
// TBD: fillChanGaps?
2441
+
2442
+
// Now apply reference phasor to all antennas
2443
+
Matrix<Float> thissol;
2444
+
for (Int iant=0;iant<shB(2);++iant) {
2445
+
thissol.reference(solB.xyPlane(iant));
2446
+
thissol-=rB;
2447
+
}
2448
+
2449
+
// Set refant, so we can put it back
2450
+
ant2B=currrefant;
2451
+
2452
+
// put back referenced solutions
2453
+
ctiter.setfparam(solB);
2454
+
ctiter.setantenna2(ant2B);
2455
+
2456
+
// Next time thru, solB is previous
2457
+
solA.reference(solB);
2458
+
flA.reference(flB);
2459
+
ant1A.reference(ant1B);
2460
+
solB.resize(); // (break references)
2461
+
flB.resize();
2462
+
ant1B.resize();
2463
+
2464
+
lastrefant=currrefant;
2465
+
first=false; // avoid first-pass stuff from now on
2466
+
2467
+
} // found
2468
+
else {
2469
+
logSink()
2470
+
<< "At "
2471
+
<< MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2472
+
<< " ("
2473
+
<< "Spw=" << ctiter.thisSpw()
2474
+
<< ", Fld=" << ctiter.thisField()
2475
+
<< ")"
2476
+
<< ", refant (id=" << currrefant
2477
+
<< ") was flagged; flagging all antennas strictly."
2478
+
<< LogIO::POST;
2479
+
// Flag all solutions in this interval
2480
+
flB.set(True);
2481
+
ctiter.setflag(flB);
2482
+
}
2483
+
2484
+
// advance to the next interval
2485
+
lastspw=ispw;
2486
+
ctiter.next();
2487
+
}
2488
+
2489
+
if (usedaltrefant)
2490
+
logSink() << LogIO::NORMAL
2491
+
<< " NB: An alternate refant was used at least once to maintain" << endl
2492
+
<< " phase continuity where the user's preferred refant drops out." << endl
2493
+
<< " Alternate refants are held constant in phase (_not_ zeroed)" << endl
2494
+
<< " during these periods, and the preferred refant may return at" << endl
2495
+
<< " a non-zero phase. This is generally harmless."
2496
+
<< LogIO::POST;
2497
+
2498
+
return;
2499
+
2500
+
}
2501
+
2188
2502
} //# NAMESPACE CASA - END
2189
2503
2190
2504
2191
2505
/*
2192
2506
Example of use in at Python level:
2193
2507
2194
2508
fringefit(vis="n14c2.ms", caltable="fail.fj", field="",spw="1",intent="",
2195
2509
selectdata=True, timerange="", antenna="", scan="5", observation="",
2196
2510
msselect="", solint="inf", refant="EF", minsnr=3.0, append=False,
2197
2511
gaintable=['n14c2.gcal'], parang=False)