Line data Source code
1 : //# GJonesPoly.cc: Implementation of GJonesPoly.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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 : //# $Id: GJonesPoly.cc,v 19.15 2006/02/03 00:29:52 gmoellen Exp $
27 :
28 : #include <synthesis/MeasurementComponents/GSpline.h>
29 : #include <synthesis/MeasurementEquations/VisEquation.h>
30 :
31 : #include <casacore/casa/Logging/LogIO.h>
32 : #include <casacore/casa/Utilities/Assert.h>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
36 : #include <casacore/casa/Quanta/MVTime.h>
37 : //#include <casa/Containers/SimOrdMap.h>
38 : #include <casacore/casa/Containers/Block.h>
39 : #include <cmath>
40 : #include <casacore/casa/fstream.h>
41 :
42 : #include <casacore/casa/System/PGPlotter.h>
43 : #include <casacore/ms/MeasurementSets/MSColumns.h>
44 : #include <msvis/MSVis/VisBuffAccumulator.h>
45 : #include <synthesis/CalTables/GJonesMBuf.h>
46 : #include <synthesis/CalTables/GJonesTable.h>
47 : #include <synthesis/CalTables/CalIter.h>
48 : #include <casacore/tables/Tables/TableUtil.h>
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : using casacore::operator*;
54 :
55 : // Define external CLIC solvers
56 : #define NEED_UNDERSCORES
57 : #if defined(NEED_UNDERSCORES)
58 : #define polyant polyant_
59 : #define splinant splinant_
60 : #define getbspl getbspl_
61 : #define phaseant phaseant_
62 : #define ampliant ampliant_
63 : #define cheb cheb_
64 : #endif
65 :
66 : extern "C" {
67 : void polyant(Int*, Int*, Int*, Int*, Int*, Int*, Int*, Int*,
68 : Double*, Double*, Double*, Double*, Double*, Double*,
69 : Double*, Double* );
70 : void splinant(Int*, Int*, Int*, Int*, Int*, Int*, Int*, Int*,
71 : Double*, Double*, Double*, Double*, Double*, Double*,
72 : Double*, Double*, Double* );
73 :
74 : void getbspl(Int*, Double*, Double*, Double*, Double*, Int*);
75 :
76 : void phaseant(Int*, Int*, Int*, Int*, Double*, Double*, Double*, Double*,
77 : Int*, Double*, Double*);
78 :
79 : void ampliant(Int*, Int*, Int*, Int*, Double*, Double*, Double*, Double*,
80 : Int*, Double*, Double*);
81 :
82 : void cheb(Int*, Double*, Double*, Int*);
83 : }
84 :
85 : //----------------------------------------------------------------------------
86 :
87 0 : GJonesSpline::GJonesSpline (VisSet& vs) :
88 : VisCal(vs),
89 : VisMueller(vs),
90 : GJones(vs),
91 0 : vs_p(&vs),
92 0 : solveAmp_p(false),
93 0 : solvePhase_p(true),
94 0 : splinetime_p(7200.0),
95 : //cacheTimeValid_p(0),
96 0 : calBuffer_p(NULL),
97 : //rawPhaseRemoval_p(false),
98 0 : solTimeStamp_p(0.0)
99 : {
100 : // Construct from a visibility set
101 : // Input:
102 : // vs VisSet& Visibility set
103 : // Output to private data:
104 : // solveAmp_p Bool true if mode_p includes amp. soln.
105 : // solvePhase_p Bool true if mode_p includes phase soln.
106 : // cacheTimeValid_p Double Time for which the current
107 : // calibration cache is valid
108 : // calBuffer_p GJonesSplineMBuf* Ptr to the applied cal. buffer
109 : //
110 :
111 0 : if (prtlev()>2) cout << "GSpline::GSpline(vs)" << endl;
112 :
113 : // Mark the Jones matrix as neither solved for nor applied,
114 : // pending initialization by setSolver() or setInterpolation()
115 0 : setSolved(false);
116 0 : setApplied(false);
117 :
118 0 : };
119 :
120 0 : GJonesSpline::GJonesSpline (const MSMetaInfoForCal& msmc) :
121 : VisCal(msmc),
122 : VisMueller(msmc),
123 : GJones(msmc),
124 0 : vs_p(NULL),
125 0 : solveAmp_p(false),
126 0 : solvePhase_p(true),
127 0 : splinetime_p(7200.0),
128 : // cacheTimeValid_p(0),
129 0 : calBuffer_p(NULL),
130 : // rawPhaseRemoval_p(false),
131 0 : solTimeStamp_p(0.0)
132 : {
133 : // Construct from a MSMetaInfoForCal
134 : // Input:
135 : // msmc MSMetaInfoForCal& MS Meta info server
136 : // Output to private data:
137 : // solveAmp_p Bool true if mode_p includes amp. soln.
138 : // solvePhase_p Bool true if mode_p includes phase soln.
139 : // cacheTimeValid_p Double Time for which the current
140 : // calibration cache is valid
141 : // calBuffer_p GJonesSplineMBuf* Ptr to the applied cal. buffer
142 : //
143 :
144 0 : if (prtlev()>2) cout << "GSpline::GSpline(msmc)" << endl;
145 :
146 : // Mark the Jones matrix as neither solved for nor applied,
147 : // pending initialization by setSolver() or setInterpolation()
148 0 : setSolved(false);
149 0 : setApplied(false);
150 :
151 0 : };
152 :
153 : //----------------------------------------------------------------------------
154 :
155 0 : GJonesSpline::~GJonesSpline ()
156 : {
157 : // Virtual destructor
158 : // Output to private data:
159 : // calBuffer_p GJonesSplineMBuf* Ptr to the applied cal. buffer
160 : //
161 :
162 0 : if (prtlev()>2) cout << "GSpline::~GSpline(vs)" << endl;
163 :
164 : // Delete the calibration buffer
165 0 : if (calBuffer_p) delete(calBuffer_p);
166 0 : };
167 :
168 : //----------------------------------------------------------------------------
169 :
170 0 : void GJonesSpline::setSolve(const Record& solvepar)
171 : {
172 : // Set the solver parameters
173 : // Input:
174 : // solvepar const Record& Solver parameters
175 : // Output to private data:
176 : // splinetime_p Double Spline knot timescale
177 :
178 0 : if (prtlev()>2) cout << "GSpline::setSolve()" << endl;
179 :
180 : // Call parent for generic pars
181 0 : SolvableVisCal::setSolve(solvepar);
182 :
183 : // Total solution interval is always all selecte data (for now)
184 0 : interval()=DBL_MAX;
185 :
186 : // Override nominal preavg handling
187 : // (avoids interpretting -1 as pre-avg up to full interval)
188 0 : if (solvepar.isDefined("preavg"))
189 0 : preavg() = solvepar.asDouble("preavg");
190 :
191 : // Spline-specific pars:
192 0 : if (solvepar.isDefined("splinetime"))
193 0 : splinetime_p = solvepar.asDouble("splinetime");
194 0 : if (solvepar.isDefined("numpoint"))
195 0 : numpoint_p = solvepar.asInt("numpoint");
196 0 : if (solvepar.isDefined("phasewrap")) // deg->rad
197 0 : phaseWrap_p = solvepar.asDouble("phasewrap")*C::pi/180.0;
198 :
199 : // Interpret mode for SPLINE case
200 : // (apmode now set in SVC::setsolve)
201 0 : solveAmp_p = (apmode().contains("A"));
202 0 : solvePhase_p = (apmode().contains("P"));
203 :
204 : // cout << "solveAmp_p = " << solveAmp_p << endl;
205 : // cout << "solvePhase_p = " << solvePhase_p << endl;
206 : // cout << "calTableName() = " << calTableName() << endl;
207 :
208 : // Mark the Jones matrix as being solved for
209 0 : setSolved(true);
210 :
211 0 : return;
212 : };
213 :
214 : //----------------------------------------------------------------------------
215 :
216 0 : void GJonesSpline::setApply(const Record& applypar)
217 : {
218 : // Set the apply parameters
219 : // Input:
220 : // applypar const Record& Interpolation parameters
221 : // Output to private data:
222 : // applyTable_p String Cal. table name containing
223 : // solutions to be applied
224 : // applySelect_p String Selection for the applied
225 : // calibration table
226 : // applyInterval_p Double Interpolation interval
227 : //
228 :
229 0 : if (prtlev()>2) cout << "GSpline::setApply()" << endl;
230 :
231 : // Extract the parameters
232 0 : if (applypar.isDefined("table"))
233 0 : calTableName() = applypar.asString("table");
234 0 : if (applypar.isDefined("select"))
235 0 : calTableSelect() = applypar.asString("select");
236 0 : if (applypar.isDefined("t"))
237 0 : interval() = applypar.asDouble("t");
238 :
239 : // Attach a calibration buffer and iterator to the calibration
240 : // table containing corrections to be applied
241 0 : GJonesSplineTable calTable(calTableName(), Table::Update);
242 0 : CalIter calIter(calTable);
243 : // Create the buffer and synchronize with the iterator
244 0 : if (calBuffer_p) delete calBuffer_p;
245 0 : calBuffer_p = new GJonesSplineMBuf(calIter);
246 0 : calBuffer_p->synchronize();
247 0 : calBuffer_p->fillCache();
248 :
249 : // Mark the Jones matrix as being applied
250 0 : setApplied(true);
251 :
252 0 : return;
253 0 : };
254 :
255 : //----------------------------------------------------------------------------
256 :
257 0 : void GJonesSpline::selfGatherAndSolve (VisSet& vs, VisEquation& ve)
258 : {
259 : // Solver for the electronic gain in spline form
260 : // Input:
261 : // me VisEquation& Measurement Equation (ME) in
262 : // which this Jones matrix resides
263 : // Output:
264 : // solve Bool true is solution succeeded
265 : // else false
266 : //
267 :
268 0 : if (prtlev()>2) cout << "GSpline::selfGatherAndSolve(vs,ve)" << endl;
269 :
270 :
271 : // cout << "Entering GSpline::solve." << endl;
272 :
273 :
274 0 : LogIO os (LogOrigin("GJonesSpline", "solve()", WHERE));
275 :
276 :
277 : os << LogIO::NORMAL
278 : << "Fitting time-dependent cubic splines."
279 0 : << LogIO::POST;
280 0 : if (solvePhase_p) {
281 : os << LogIO::NORMAL
282 : << "Solving for phase splines with splinetime= " << splinetime_p
283 0 : << LogIO::POST;
284 : }
285 0 : if (solveAmp_p) {
286 : os << LogIO::NORMAL
287 : << "Solving for amplitude splines with splinetime= " << splinetime_p
288 0 : << LogIO::POST;
289 : }
290 :
291 : // Arrange for iteration over data
292 0 : Block<Int> columns;
293 : // avoid scan iteration
294 0 : columns.resize(4);
295 0 : columns[0]=MS::ARRAY_ID;
296 0 : columns[1]=MS::FIELD_ID;
297 0 : columns[2]=MS::DATA_DESC_ID;
298 0 : columns[3]=MS::TIME;
299 0 : vs.resetVisIter(columns,interval());
300 0 : VisIter& vi(vs.iter());
301 0 : VisBuffer vb(vi);
302 :
303 : // Count polarizations
304 0 : Int nPH(0);
305 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
306 0 : vi.origin();
307 0 : Int ncorr(vb.corrType().nelements());
308 0 : nPH= max(nPH,min(ncorr,2));
309 : }
310 :
311 : // nPH has the number of correlations (1 or 2) to process
312 :
313 : // Initialize time-series accumulation buffers for the
314 : // corrected and corrupted visibility data and associated
315 : // weights.
316 0 : std::map<String,Int> timeValueMap;
317 0 : Vector<Double> timeValues;
318 0 : PtrBlock<Matrix<Complex>* > visTimeSeries;
319 0 : PtrBlock<Matrix<Double>* > weightTimeSeries;
320 0 : Int nTimeSeries = 0;
321 :
322 : // Initialize antenna indices
323 0 : Vector<Int> ant1(nBln(), -1);
324 0 : Vector<Int> ant2(nBln(), -1);
325 :
326 0 : Int k=0;
327 0 : for (Int a1=0; a1 < nAnt(); a1++) {
328 0 : for (Int a2=a1; a2 < nAnt(); a2++,k++) {
329 0 : ant1(k) = a1 + 1;
330 0 : ant2(k) = a2 + 1;
331 : // cout << a1 << " " << a2 << " " << k << " " << blnidx(a1,a2) << endl;
332 : };
333 : };
334 :
335 : // Iterate, accumulating the averaged visibility time-series
336 : Int chunk;
337 : // Mean reference frequency
338 0 : Double meanFreq = 0;
339 0 : Int nMeanFreq = 0;
340 :
341 : // With current VisIter sort order, this chunking does
342 : // all times for each spw and field
343 : // Double t0;
344 0 : for (chunk=0, vi.originChunks(); vi.moreChunks(); vi.nextChunk(), chunk++) {
345 :
346 : // Extract the current visibility buffer spectral window id.
347 : // and number for frequency channels
348 0 : Int spwid = vi.spectralWindow();
349 0 : Int nChan = vs.numberChan()(spwid);
350 0 : String fieldName = vi.fieldName();
351 :
352 : os << LogIO::NORMAL
353 : << "Accumulating data for: field= " << fieldName
354 : << ", spw= " << spwid
355 : << ", nchan= " << nChan
356 0 : << LogIO::POST;
357 :
358 : // Compute the corrected and corrupted data at the position of
359 : // this Jones matrix in the Measurement Equation. The corrupted
360 : // data are the model visibilities propagated along the ME from
361 : // the sky to the immediate right of the current Jones matrix.
362 : // The corrected data are the observed data corrected for all
363 : // Jones matrices up to the immediate left of the current Jones
364 : // matrix.
365 :
366 : // Arrange to accumulate
367 0 : VisBuffAccumulator vba(nAnt(),preavg(),false);
368 :
369 0 : vi.origin();
370 : // t0=86400.0*floor(vb.time()(0)/86400.0);
371 :
372 : // Collapse each timestamp in this chunk according to VisEq
373 : // with calibration and averaging
374 0 : for (vi.origin(); vi.more(); vi++) {
375 :
376 : // TBD: initialize weights from sigma here?
377 :
378 0 : ve.collapse(vb);
379 :
380 : // If permitted/required by solvable component, normalize
381 0 : if (normalizable())
382 0 : vb.normalize();
383 :
384 : // If this solve not freqdep, and channels not averaged yet, do so
385 0 : if (!freqDepMat() && vb.nChannel()>1)
386 0 : vb.freqAveCubes();
387 :
388 : // Accumulate collapsed vb in a time average
389 0 : vba.accumulate(vb);
390 : }
391 0 : vba.finalizeAverage();
392 :
393 : // The VisBuffer to work with in solve
394 0 : VisBuffer& svb(vba.aveVisBuff());
395 :
396 : // cout << " (Chunk accumulated) nrow = " << svb.nRow() << endl;
397 : // cout << " nChan = " << svb.nChannel() << endl;
398 :
399 : // NB: The svb VisBuffer has many timestamps in it....
400 :
401 : // index to parallel hands (assumes canonical order)
402 0 : Int nCorr = svb.corrType().nelements();
403 0 : Vector<Int> polidx(2,0);
404 0 : polidx(1)=nCorr-1;
405 :
406 : // Accumulate mean frequency of the averaged data
407 0 : meanFreq += mean(svb.frequency());
408 0 : nMeanFreq++;
409 :
410 : // Iterate over the current accumulated visibility buffer, obtaining
411 : // a common time series in the visibility data
412 0 : for (Int row=0; row < svb.nRow(); row++) {
413 : // Antenna numbers
414 0 : Int ant1num = svb.antenna1()(row);
415 0 : Int ant2num = svb.antenna2()(row);
416 :
417 : // Reject auto-correlation data
418 0 : if (ant1num != ant2num) {
419 : // Compute baseline index
420 0 : Int baselineIndex = blnidx(ant1num,ant2num);
421 :
422 : // Weight
423 0 : Vector<Float> rowWeight(svb.weightMat().column(row));
424 0 : if (sum(rowWeight) > 0) {
425 : // Map the current time stamp to a time series index
426 : // using 0.1 second precision, i.e. time values within
427 : // a tenth of a second of each other will be accumulated
428 : // within the same bin.
429 0 : MVTime mvt(svb.time()(row) / C::day);
430 0 : String timeKey = mvt.string(MVTime::TIME, 7);
431 0 : Int timeIndex = 0;
432 : // Check the time stamp index to this precision
433 0 : if (timeValueMap.find(timeKey) != timeValueMap.end( )) {
434 0 : timeIndex = timeValueMap[timeKey];
435 : } else {
436 : // Create a new time series entry
437 0 : timeIndex = nTimeSeries++;
438 0 : timeValueMap.insert(std::pair<casacore::String, casacore::Int>(timeKey, timeIndex));
439 0 : timeValues.resize(nTimeSeries, true);
440 0 : timeValues(timeIndex) = svb.time()(row);
441 0 : Complex czero(0,0);
442 0 : visTimeSeries.resize(nTimeSeries, true);
443 0 : visTimeSeries[timeIndex] = new Matrix<Complex>(nBln(),nPH, czero);
444 0 : weightTimeSeries.resize(nTimeSeries, true);
445 0 : weightTimeSeries[timeIndex] = new Matrix<Double>(nBln(),nPH, 0);
446 : };
447 :
448 : // Accumulate the current visbility row in the common time series
449 0 : for (Int ip=0;ip<nPH;++ip) {
450 0 : Double dwt=Double(rowWeight(polidx(ip)));
451 0 : (*visTimeSeries[timeIndex])(baselineIndex,ip) +=
452 0 : svb.visCube()(polidx(ip),0,row) * dwt;
453 0 : (*weightTimeSeries[timeIndex])(baselineIndex,ip) += dwt;
454 : }
455 0 : };
456 0 : };
457 : };
458 0 : }; // for (chunk...) iteration
459 :
460 : // Create amplitude, phase and weight arrays per time-slot
461 : // and interferometer index in the form required by the
462 : // GILDAS solver, splinant.
463 0 : Int nTimes = timeValueMap.size( );
464 :
465 : os << LogIO::NORMAL
466 : << "Number of timestamps in data = " << nTimes
467 0 : << LogIO::POST;
468 :
469 0 : Vector<Double> time(nTimes, 0);
470 0 : Cube<Double> amp(nTimes, nBln(), nPH, 0);
471 0 : Cube<Double> phase(nTimes, nBln(), nPH, 0);
472 0 : Cube<Double> weight(nTimes, nBln(),nPH, 0);
473 : //
474 : // First create a simple index in order of ascending time ordinate
475 0 : Vector<Int> index(nTimes);
476 0 : indgen(index);
477 0 : for (Int j=0; j < nTimes; j++) {
478 0 : for (Int k=j+1; k < nTimes; k++) {
479 0 : if (timeValues(index(j)) > timeValues(index(k))) {
480 0 : Int swap = index(j);
481 0 : index(j) = index(k);
482 0 : index(k) = swap;
483 : };
484 : };
485 : };
486 :
487 : // Fill the amplitude, phase and weight arrays with normalized
488 : // and scaled vaues from the accumulated time series.
489 0 : for (Int t=0; t < nTimes; t++) {
490 0 : for (Int baselineIndex=0; baselineIndex < nBln(); ++baselineIndex) {
491 : // Iterate in ascending time order
492 0 : Int indx = index(t);
493 : // Time
494 0 : time(t) = timeValues(indx);
495 :
496 0 : for (Int ip=0;ip<nPH;++ip) {
497 :
498 0 : Double wgt = (*weightTimeSeries[indx])(baselineIndex,ip);
499 0 : if (wgt > 0) {
500 0 : Complex vis = (*visTimeSeries[indx])(baselineIndex,ip);
501 :
502 : // Amplitude
503 0 : if (abs(vis) > 0) {
504 0 : amp(t,baselineIndex,ip) = log(abs(vis)/wgt);
505 : };
506 :
507 : // Phase
508 0 : if(solvePhase_p){
509 0 : Double visPhase = arg(vis);
510 0 : visPhase = remainder(visPhase, 2*C::pi);
511 0 : if (visPhase > C::pi) {
512 0 : visPhase -= 2 * C::pi;
513 0 : } else if (visPhase < -C::pi) {
514 0 : visPhase += 2 * C::pi;
515 : };
516 0 : phase(t,baselineIndex,ip) = visPhase;
517 : }
518 :
519 :
520 : // Weight
521 0 : weight(t,baselineIndex,ip) = wgt;
522 : //weight(t,baselineIndex,ip) = 1.0;
523 : }
524 : };
525 : };
526 : };
527 :
528 : // cout << "weight = " << weight << endl;
529 : // cout << "weight.shape() = " << weight.shape() << endl;
530 :
531 : // Delete the accumulation buffers
532 0 : for (Int k=0; k < nTimes; k++) {
533 0 : delete(visTimeSeries[k]);
534 0 : delete(weightTimeSeries[k]);
535 : };
536 :
537 :
538 : // Double t0=86400.0*floor(min(time)/86400.0);
539 :
540 : // Compute the number of spline knots
541 0 : Vector<Double> knots, ampKnots(1,0.0), phaKnots(1,0.0);
542 0 : Int nKnots = getKnots(time, knots);
543 :
544 : os << LogIO::NORMAL
545 : << "Number of cubic spline control points = " << nKnots-4
546 0 : << LogIO::POST;
547 :
548 : os << LogIO::NORMAL
549 : << "Number of cubic spline knots = " << nKnots
550 0 : << LogIO::POST;
551 :
552 : os << LogIO::NORMAL
553 : << "Number of cubic spline segments = " << nKnots-7
554 0 : << LogIO::POST;
555 :
556 0 : if (nTimes < (nKnots-4) ) {
557 : os << LogIO::SEVERE
558 : << "Insufficient number of timestamps for cubic splines:" << endl
559 0 : << " Control points - nTimes > 0 (" << nKnots-4-nTimes << ")"
560 0 : << LogIO::POST;
561 : os << LogIO::SEVERE
562 : << "Increase splinetime or reduce preavg (if possible), or " << endl
563 : << " solve for ordinary G instead."
564 0 : << LogIO::POST;
565 0 : return;
566 : };
567 :
568 0 : solTimeStamp_p = mean(knots);
569 :
570 : // Declare work arrays and returned value
571 : // arrays to be used by the solver
572 0 : Matrix<Double> wk1(4, nTimes);
573 0 : Matrix<Double> wk2(4*nAnt(), nKnots*nAnt());
574 0 : Vector<Double> wk3(nKnots*nAnt());
575 0 : Vector<Double> rmsfit(nBln(),0);
576 :
577 0 : Double dzero(0);
578 :
579 : // GSpline is nominally dual-pol
580 0 : Cube<Double> polyCoeffAmp(nAnt(),nKnots,2, dzero);
581 0 : Cube<Double> polyCoeffPhase(nAnt(),nKnots,2, dzero);
582 :
583 : // Fit using spline polynomials (GILDAS splinant)
584 : Bool dum;
585 : Int iy;
586 : // GILDAS solver uses one-relative antenna numbers
587 0 : Int rAnt = refant() + 1;
588 0 : Int nBL(nBln());
589 0 : Int nA(nAnt());
590 :
591 0 : if (solveAmp_p) {
592 : // Amplitude solution
593 :
594 : os << LogIO::NORMAL
595 : << "Fitting amplitude spline."
596 0 : << LogIO::POST;
597 :
598 0 : ampKnots.resize(knots.shape());
599 0 : ampKnots=knots;
600 :
601 0 : iy = 1;
602 :
603 0 : for (Int ip=0;ip<nPH;++ip) {
604 :
605 0 : splinant(&iy,
606 : &nTimes,
607 : &nBL,
608 : ant1.getStorage(dum),
609 : ant2.getStorage(dum),
610 : &rAnt,
611 : &nKnots,
612 : &nA,
613 : time.getStorage(dum),
614 0 : amp.xyPlane(ip).getStorage(dum),
615 0 : weight.xyPlane(ip).getStorage(dum),
616 : knots.getStorage(dum),
617 : wk1.getStorage(dum),
618 : wk2.getStorage(dum),
619 : wk3.getStorage(dum),
620 : rmsfit.getStorage(dum),
621 0 : polyCoeffAmp.xyPlane(ip).getStorage(dum));
622 :
623 0 : ostringstream o;
624 0 : o << calTableName() << ".AMP.pol" << ip << ".log";
625 0 : String logfile=o.str();
626 0 : writeAsciiLog(logfile, polyCoeffAmp.xyPlane(ip), rmsfit, false);
627 : // plotsolve(time, amp, weight, rmsfit, polyCoeffAmp, false);
628 :
629 0 : }
630 :
631 : }
632 :
633 0 : if (solvePhase_p){
634 : // Phase solution
635 :
636 : os << LogIO::NORMAL
637 : << "Searching for and correcting phase-wraps on each baseline."
638 0 : << LogIO::POST;
639 :
640 0 : for (Int ip=0;ip<nPH;++ip) {
641 :
642 0 : for(Int bsInd=0; bsInd < nBln() ; ++bsInd){
643 :
644 0 : Int npoi=numpoint_p;
645 0 : if (npoi < 2) npoi=2;
646 :
647 0 : for (Int k1=npoi; k1 < nTimes-npoi; k1=k1+(npoi/2)){
648 0 : Vector<Double> vecAheadPh(npoi), vecBheadPh(npoi);
649 0 : for (Int k=0; k < npoi; ++k){
650 0 : vecBheadPh[k]=phase(k1-k-1,bsInd,ip);
651 0 : vecAheadPh[k]=phase(k1+k, bsInd,ip);
652 : }
653 0 : Double medA=median(vecAheadPh);
654 0 : Double medB=median(vecBheadPh);
655 0 : if((medA- medB) >phaseWrap_p)
656 0 : for(Int k=k1; k< nTimes; ++k)
657 0 : phase(k, bsInd,ip) -= 2*C::pi;
658 0 : if((medA-medB) < -phaseWrap_p)
659 0 : for(Int k=k1; k< nTimes; ++k)
660 0 : phase(k, bsInd,ip) += 2*C::pi;
661 :
662 0 : }
663 :
664 0 : if(nTimes >npoi*2){
665 :
666 : //the last npoi pts
667 0 : Vector<Double> vecAheadPh(npoi), vecBheadPh(npoi);
668 0 : for (Int k=0; k <npoi; ++k){
669 0 : vecBheadPh[k]=phase(nTimes-2*npoi+k,bsInd,ip);
670 0 : vecAheadPh[k]=phase(nTimes-npoi+k, bsInd,ip);
671 : }
672 0 : Double medA=median(vecAheadPh);
673 0 : Double medB=median(vecBheadPh);
674 0 : if((medA- medB) > phaseWrap_p)
675 0 : for(Int k=(nTimes-npoi); k< nTimes; ++k)
676 0 : phase(k, bsInd,ip) -= 2*C::pi;
677 0 : if((medA-medB) < -phaseWrap_p)
678 0 : for(Int k=(nTimes-npoi); k< nTimes; ++k)
679 0 : phase(k, bsInd,ip) += 2*C::pi;
680 :
681 : // The first npoi points
682 :
683 0 : for (Int k=0; k <npoi; ++k){
684 0 : vecBheadPh[k]=phase(k,bsInd,ip);
685 0 : vecAheadPh[k]=phase(npoi+k, bsInd,ip);
686 : }
687 0 : medA=median(vecAheadPh);
688 0 : medB=median(vecBheadPh);
689 0 : if((medA- medB) > phaseWrap_p)
690 0 : for(Int k=0; k< npoi; ++k)
691 0 : phase(k, bsInd,ip) += 2*C::pi;
692 0 : if((medA-medB) < -phaseWrap_p)
693 0 : for(Int k=0; k< npoi; ++k)
694 0 : phase(k, bsInd,ip) -= 2*C::pi;
695 :
696 0 : }
697 : }
698 : }
699 :
700 : os << LogIO::NORMAL
701 : << "Fitting phase spline."
702 0 : << LogIO::POST;
703 :
704 0 : phaKnots.resize(knots.shape());
705 0 : phaKnots=knots;
706 0 : iy = 2;
707 :
708 0 : for (Int ip=0;ip<nPH;++ip) {
709 :
710 0 : wk1=0.0; wk2=0.0; wk3=0.0;
711 :
712 0 : splinant(&iy,
713 : &nTimes,
714 : &nBL,
715 : ant1.getStorage(dum),
716 : ant2.getStorage(dum),
717 : &rAnt,
718 : &nKnots,
719 : &nA,
720 : time.getStorage(dum),
721 0 : phase.xyPlane(ip).getStorage(dum),
722 0 : weight.xyPlane(ip).getStorage(dum),
723 : knots.getStorage(dum),
724 : wk1.getStorage(dum),
725 : wk2.getStorage(dum),
726 : wk3.getStorage(dum),
727 : rmsfit.getStorage(dum),
728 0 : polyCoeffPhase.xyPlane(ip).getStorage(dum));
729 :
730 :
731 : // cout << "Finished solve." << endl;
732 :
733 :
734 0 : ostringstream o;
735 0 : o << calTableName() << ".PHASE.pol" << ip << ".log";
736 0 : String logfile=o.str();
737 0 : writeAsciiLog(logfile, polyCoeffPhase.xyPlane(ip), rmsfit, true);
738 :
739 : // TBD: make multi-pol plotsolve
740 : // plotsolve(time, phase, weight, rmsfit, polyCoeffPhase, true);
741 :
742 0 : }
743 :
744 :
745 :
746 : };
747 :
748 : // Update the output calibration table
749 0 : Vector<Int> antId(nAnt());
750 0 : indgen(antId);
751 0 : Vector<Int> fieldIdKeys = fieldIdRange();
752 0 : Vector<String> freqGrpName(nAnt(), "");
753 0 : Vector<String> polyType(nAnt(), "SPLINE");
754 0 : Vector<String> polyMode(nAnt(), apmode());
755 0 : if (apmode()=="AP")
756 0 : polyMode="PHASAMP";
757 0 : if (apmode()=="P")
758 0 : polyMode="PHAS";
759 0 : if (apmode()=="A")
760 0 : polyMode="AMP";
761 :
762 0 : Vector<Complex> scaleFactor(nAnt(), Complex(1,0));
763 0 : Vector<String> phaseUnits(nAnt(), "RADIANS");
764 0 : Vector<Int> refAnt(nAnt(), refant());
765 :
766 : // Use mean frequency as the reference frequency
767 0 : if (nMeanFreq > 0) meanFreq = meanFreq / nMeanFreq;
768 0 : Vector<MFrequency> refFreq(nAnt(), MFrequency(Quantity(meanFreq, "Hz")));
769 :
770 : // If not incremental calibration then create an empty
771 : // calibration buffer to hold the output calibration solutions
772 0 : if (!calBuffer_p) {
773 0 : newCalBuffer(fieldIdKeys, antId);
774 : } else {
775 : // Force an explicit read to cache for all calibration columns
776 : // (in case the input calibration table is being overwritten
777 : // in place)
778 0 : calBuffer_p->fillCache();
779 :
780 : };
781 :
782 : // Update the calibration table
783 : // cout << "ampKnots.nelements() = " << ampKnots.nelements() << endl;
784 : // cout << "phaKnots.nelements() = " << phaKnots.nelements() << endl;
785 : // cout << "nKnots = " << nKnots << endl;
786 :
787 : // cout << "polyCoeffPhase.shape() = " << polyCoeffPhase.shape() << endl;
788 : // cout << "polyCoeffPhase.reform(IPosition(nAnt(),ampKnots.nelements()*2)).shape() = "
789 : // << polyCoeffPhase.reform(IPosition(2,nAnt(),nKnots*2)).shape() << endl;
790 :
791 :
792 0 : Matrix<Double> ampco(polyCoeffAmp.reform(IPosition(2,nAnt(),nKnots*2)));
793 0 : Matrix<Double> phaco(polyCoeffPhase.reform(IPosition(2,nAnt(),nKnots*2)));
794 :
795 0 : updateCalTable(fieldIdKeys, antId, freqGrpName, polyType, polyMode,
796 : scaleFactor,
797 : ampco,phaco,
798 : // polyCoeffAmp.xyPlane(0), polyCoeffPhase.xyPlane(0),
799 : phaseUnits, ampKnots, phaKnots, refFreq, refAnt);
800 :
801 0 : return;
802 0 : };
803 :
804 : //----------------------------------------------------------------------------
805 :
806 0 : void GJonesSpline::calcPar() {
807 :
808 : // Calculate new gain parameters (in this case, the G matrix elements)
809 :
810 0 : if (prtlev()>6) cout << " GSpline::calcPar()" << endl;
811 :
812 0 : LogIO os (LogOrigin("GJonesSpline", "calcPar", WHERE));
813 :
814 0 : Double freqRatio=1.0;
815 0 : Double vbFreqHz = 1.0e9*mean(currFreq());
816 :
817 : // cout << "vbFreqHz = " << vbFreqHz << endl;
818 :
819 0 : currCPar().resize(nPar(),1,nAnt());
820 0 : currParOK().resize(nPar(),1,nAnt());
821 0 : currParOK()=false;
822 :
823 : // Compute them, per antenna
824 0 : for (Int iant=0; iant < nAnt(); iant++) {
825 :
826 : // Match this antenna id. and the visibility buffer field id.
827 : // in the calibration buffer
828 : Vector<Int> matchingRows =
829 0 : calBuffer_p->matchAntenna1AndFieldId(iant,currField());
830 :
831 0 : if (matchingRows.nelements() > 0) {
832 : // Matching calibration solution found
833 0 : Int row = matchingRows(0);
834 0 : Vector<Double> ampVal(2,1.0);
835 0 : Vector<Double> phaseVal(2,0.0);
836 0 : Complex gain(calBuffer_p->scaleFactor()(row));
837 0 : String mode = calBuffer_p->polyMode()(row);
838 :
839 : // cout << "gain = " << gain << endl;
840 :
841 : // Compute the ratio between the calibration solution
842 : // reference frequency and the mean observed frequency
843 : // of the visibility data to be corrected
844 0 : IPosition refFreqPos = calBuffer_p->refFreqMeas().shape();
845 0 : refFreqPos = 0;
846 0 : refFreqPos.setLast(IPosition(1,row));
847 0 : MFrequency refFreq = calBuffer_p->refFreqMeas()(refFreqPos);
848 0 : Double refFreqHz = refFreq.get("Hz").getValue();
849 0 : freqRatio = abs(refFreqHz) > 0 ? vbFreqHz / refFreqHz : 1.0;
850 : // cout << "freqRatio = " << freqRatio << endl;
851 :
852 : // Compute amplitude polynomial
853 0 : if (mode.contains("AMP") || mode.contains("A&P")) {
854 : // Extract amplitude spline polynomial knots
855 0 : IPosition ampKnotsPos = calBuffer_p->splineKnotsAmp().shape();
856 0 : ampKnotsPos = 0;
857 0 : Int nAmpKnotsPos = ampKnotsPos.nelements();
858 0 : ampKnotsPos[nAmpKnotsPos-1] = row;
859 0 : Int nKnots = calBuffer_p->nKnotsAmp()(row);
860 0 : Vector<Double> ampKnots(nKnots);
861 0 : for (Int k=0; k < nKnots; k++) {
862 0 : ampKnotsPos[nAmpKnotsPos-2] = k;
863 0 : ampKnots(k) = calBuffer_p->splineKnotsAmp()(ampKnotsPos);
864 : };
865 :
866 : // Extract amplitude spline polynomial coefficients
867 0 : IPosition ampCoeffPos = calBuffer_p->polyCoeffAmp().shape();
868 0 : ampCoeffPos = 0;
869 0 : Int nAmpCoeffPos = ampCoeffPos.nelements();
870 0 : ampCoeffPos[nAmpCoeffPos-1] = row;
871 0 : Int nPoly = calBuffer_p->nPolyAmp()(row);
872 0 : Vector<Double> ampCoeff(2*nPoly);
873 0 : for (Int k=0; k < 2*nPoly; k++) {
874 0 : ampCoeffPos[nAmpCoeffPos-2] = k;
875 0 : ampCoeff(k) = calBuffer_p->polyCoeffAmp()(ampCoeffPos);
876 : };
877 :
878 : // cout << "ampCoeff = " << ampCoeff << endl;
879 :
880 : // Compute amplitude spline polynomial value
881 0 : Vector<Double> ac;
882 0 : for (Int i=0;i<2;++i) {
883 0 : ac.reference(ampCoeff(IPosition(1,i*nPoly),
884 0 : IPosition(1,(i+1)*nPoly-1)));
885 0 : ampVal(i) *= exp(getSplineVal(currTime(), ampKnots, ac));
886 : }
887 : // cout << iant << " ampVal = " << ampVal << endl;
888 0 : };
889 :
890 : // Compute phase polynomial
891 0 : if (mode.contains("PHAS") || mode.contains("A&P")) {
892 :
893 : // Extract phase spline polynomial knots
894 0 : IPosition phaseKnotsPos = calBuffer_p->splineKnotsPhase().shape();
895 0 : phaseKnotsPos = 0;
896 0 : Int nPhaseKnotsPos = phaseKnotsPos.nelements();
897 0 : phaseKnotsPos[nPhaseKnotsPos-1] = row;
898 0 : Int nKnots = calBuffer_p->nKnotsPhase()(row);
899 0 : Vector<Double> phaseKnots(nKnots);
900 0 : for (Int k=0; k < nKnots; k++) {
901 0 : phaseKnotsPos[nPhaseKnotsPos-2] = k;
902 0 : phaseKnots(k) = calBuffer_p->splineKnotsPhase()(phaseKnotsPos);
903 : };
904 :
905 : // Extract phase spline polynomial coefficients
906 0 : IPosition phaseCoeffPos = calBuffer_p->polyCoeffPhase().shape();
907 0 : phaseCoeffPos = 0;
908 0 : Int nPhaseCoeffPos = phaseCoeffPos.nelements();
909 0 : phaseCoeffPos[nPhaseCoeffPos-1] = row;
910 0 : Int nPoly = calBuffer_p->nPolyPhase()(row);
911 0 : Vector<Double> phaseCoeff(2*nPoly);
912 0 : for (Int k=0; k < 2*nPoly; k++) {
913 0 : phaseCoeffPos[nPhaseCoeffPos-2] = k;
914 0 : phaseCoeff(k) = calBuffer_p->polyCoeffPhase()(phaseCoeffPos);
915 : };
916 :
917 : // cout << "phaseCoeff = " << phaseCoeff << endl;
918 :
919 : // Compute phase spline polynomial value
920 0 : Vector<Double> pc;
921 0 : for (Int i=0;i<2;++i) {
922 0 : pc.reference(phaseCoeff(IPosition(1,i*nPoly),
923 0 : IPosition(1,(i+1)*nPoly-1)));
924 :
925 : // cout << "pc = " << pc << endl;
926 :
927 0 : phaseVal(i) = getSplineVal(currTime(), phaseKnots, pc);
928 :
929 : // Handle gildas sign convention on spline phases
930 0 : phaseVal(i) = -phaseVal(i);
931 :
932 : // Scale by the ratio of the observing frequency of the
933 : // data to be corrected and the reference frequency of the
934 : // calibration solution
935 0 : phaseVal(i) *= freqRatio;
936 : }
937 0 : };
938 :
939 : // cout << "phaseVal = " << phaseVal << endl;
940 :
941 : // Fill the (matrix element) parameters
942 0 : for (Int i=0;i<2;++i) {
943 0 : currCPar()(i,0,iant) = gain * ampVal(i) * Complex(cos(phaseVal(i)),
944 0 : sin(phaseVal(i)) );
945 0 : currParOK()(i,0,iant) = true;
946 : }
947 :
948 0 : };
949 0 : };
950 :
951 : //cout << "currCPar() = " << currCPar() << endl;
952 :
953 0 : validateP();
954 0 : invalidateJ();
955 :
956 0 : return;
957 0 : };
958 :
959 : //----------------------------------------------------------------------------
960 :
961 0 : void GJonesSpline::newCalBuffer (const Vector<Int>& fieldIdKeys,
962 : const Vector<Int>& antennaId)
963 : {
964 : // Create and fill and empty output calibration buffer
965 : // Input:
966 : // fieldIdKeys const Vec<Int>& Field id.'s to span
967 : // antennaId const Vec<Int>& Antenna id.'s to span
968 : // Output to private data:
969 : // calBuffer_p GJonesSplineMBuf* Calibration buffer
970 : //
971 0 : LogIO os (LogOrigin("GJonesSpline", "newCalBuffer", WHERE));
972 :
973 : // Delete calibration buffer if it already exists
974 0 : if (calBuffer_p) delete calBuffer_p;
975 :
976 : // Initialize the GJonesSpline calibration buffer, spanning the
977 : // antenna id.'s and field id.'s specified.
978 0 : Vector<Int> key(2);
979 0 : key(0) = MSC::ANTENNA1;
980 0 : key(1) = MSC::FIELD_ID;
981 0 : Block<Vector<Int> > keyvals(2);
982 0 : keyvals[0] = antennaId;
983 0 : keyvals[1] = fieldIdKeys;
984 0 : calBuffer_p = new GJonesSplineMBuf(key, keyvals);
985 :
986 0 : return;
987 0 : };
988 :
989 : //----------------------------------------------------------------------------
990 :
991 0 : Int GJonesSpline::getKnots (const Vector<Double>& times, Vector<Double>& knots)
992 : {
993 : // Compute the number and location of the spline knots
994 : // Input:
995 : // times const Vector<Double>& Vector of time values
996 : // Output:
997 : // knots Vector<Double>& Knot positions
998 : // getKnots Int Number of knots
999 : //
1000 0 : LogIO os (LogOrigin("GJonesSpline", "getKnots()", WHERE));
1001 :
1002 : // Use algorithm in GILDAS, with user-defined timescale
1003 0 : Int n=times.nelements();
1004 0 : Int nSeg = Int(0.5 + (times(n-1)-times(0))/splinetime_p);
1005 0 : nSeg = max(1,nSeg); // never fewer than 1, obviously
1006 0 : Int ncenterKnots = nSeg -1;
1007 0 : Double step = (times(n-1) - times(0)) / nSeg;
1008 :
1009 : os << LogIO::NORMAL
1010 : << "Gridded splinetime = " << step << " sec."
1011 0 : << LogIO::POST;
1012 :
1013 0 : Int numOfknots = ncenterKnots + 8;
1014 : //cout << "numOfknots" << numOfknots << endl;
1015 0 : knots.resize(numOfknots);
1016 0 : knots(0)=knots(1) = knots(2) = knots(3) = times(0);
1017 0 : knots(numOfknots-1) = knots(numOfknots-2) = knots(numOfknots-3) =
1018 0 : knots(numOfknots-4)=times(n-1);
1019 0 : for(Int k=0; k < ncenterKnots; k++){
1020 0 : knots(k+4) = times(0)+ (k + 1) * step;
1021 : }
1022 :
1023 0 : return numOfknots;
1024 0 : };
1025 :
1026 : //----------------------------------------------------------------------------
1027 :
1028 0 : void GJonesSpline::updateCalTable (const Vector<Int>& fieldIdKeys,
1029 : const Vector<Int>& antennaId,
1030 : const Vector<String>& freqGrpName,
1031 : const Vector<String>& polyType,
1032 : const Vector<String>& polyMode,
1033 : const Vector<Complex>& scaleFactor,
1034 : const Matrix<Double>& polyCoeffAmp,
1035 : const Matrix<Double>& polyCoeffPhase,
1036 : const Vector<String>& phaseUnits,
1037 : const Vector<Double>& splineKnotsAmp,
1038 : const Vector<Double>& splineKnotsPhase,
1039 : const Vector<MFrequency>& refFreq,
1040 : const Vector<Int>& refAnt)
1041 : {
1042 : // Update the output calibration table
1043 : // Input:
1044 : // fieldIdKeys const Vec<Int>& Field id.'s to span
1045 : // antennaId const Vec<Int>& Antenna id. for each solution
1046 : // freqGrpName const Vec<String>& Freq. group name (per soln.)
1047 : // polyType const Vec<String>& Polynomial type (per soln.)
1048 : // polyMode const Vec<String>& Polynomial mode (per soln.)
1049 : // scaleFactor const Vec<Complex>& Polynomial scale factor
1050 : // (per solution)
1051 : // polyCoeffAmp const Matrix<Double>& Polynomial amplitude
1052 : // coefficients (per solution)
1053 : // polyCoeffPhase const Matrix<Double>& Polynomial phase coefficients
1054 : // (per solution)
1055 : // phaseUnits const Vec<String>& Phase units
1056 : // splineKnotsAmp const Vector<Double>& Amp. spline knot positions
1057 : // (same for all solutions)
1058 : // splineKnotsPhase const Vector<Double>& Phase spline knot positions
1059 : // (same for all solutions)
1060 : // refFreq const Vec<MFreq>& Reference freq. (per soln.)
1061 : // refAnt const Vec<Int>& Reference antenna (per soln.)
1062 : //
1063 0 : LogIO os (LogOrigin("GJonesSpline", "updateCalTable", WHERE));
1064 :
1065 : // Add each solution to the calibration buffer
1066 0 : for (uInt k=0; k < antennaId.nelements(); k++) {
1067 0 : for (uInt j=0; j < fieldIdKeys.nelements(); j++) {
1068 : // Find matching rows for this antenna and field id.
1069 : Vector<Int> matchingRows =
1070 0 : calBuffer_p->matchAntenna1AndFieldId(antennaId(k), fieldIdKeys(j));
1071 :
1072 : // Update the matching rows in the calibration buffer
1073 0 : calBuffer_p->fillMatchingRows(matchingRows, freqGrpName(k), polyType(k),
1074 0 : polyMode(k), scaleFactor(k),
1075 0 : polyCoeffAmp.row(k).nelements()/2,
1076 0 : polyCoeffPhase.row(k).nelements()/2,
1077 0 : polyCoeffAmp.row(k), polyCoeffPhase.row(k),
1078 0 : phaseUnits(k), splineKnotsAmp.nelements(),
1079 0 : splineKnotsPhase.nelements(),
1080 : splineKnotsAmp, splineKnotsPhase,
1081 0 : refFreq(k), refAnt(k));
1082 0 : };
1083 : };
1084 :
1085 : // calBuffer_p->fieldId().set(-1);
1086 0 : calBuffer_p->calDescId().set(0);
1087 : // cout << "Time = " << MVTime(solTimeStamp_p/C::day).string(MVTime::YMD,7) << endl;
1088 0 : calBuffer_p->timeMeas().set(MEpoch(MVEpoch(solTimeStamp_p/86400.0)));
1089 :
1090 : // Delete the output calibration table if it already exists
1091 0 : if (calTableName()!="" && TableUtil::canDeleteTable(calTableName())) {
1092 0 : TableUtil::deleteTable(calTableName());
1093 : };
1094 :
1095 : os << LogIO::NORMAL
1096 0 : << "Storing solutions in table " << calTableName()
1097 0 : << LogIO::POST;
1098 :
1099 : // Write the calibration buffer to the output calibration table
1100 0 : Table::TableOption accessMode = Table::New;
1101 0 : if (Table::isWritable(calTableName())) accessMode = Table::Update;
1102 0 : GJonesSplineTable calTable(calTableName(), accessMode);
1103 0 : calBuffer_p->append(calTable);
1104 :
1105 :
1106 : // Add CAL_DESC
1107 0 : CalDescRecord cdr;
1108 0 : cdr.defineSpwId(Vector<Int>(1,-1));
1109 0 : cdr.defineMSName(Path(msName()).baseName());
1110 0 : calTable.putRowDesc(0,cdr);
1111 :
1112 0 : return;
1113 0 : };
1114 :
1115 : //----------------------------------------------------------------------------
1116 :
1117 0 : Double GJonesSpline::getSplineVal (Double x,
1118 : Vector<Double>& knots,
1119 : Vector<Double>& coeff)
1120 : {
1121 : // Compute a polynomial spline value using the GILDAS routine getbspl
1122 : // Input:
1123 : // x Double Value at which to compute spline
1124 : // knots Vector<Double>& Knot locations
1125 : // coeff Vector<Double>& Spline coefficients
1126 : // Output:
1127 : // getSplineVal Double Computed spline value
1128 : //
1129 0 : LogIO os (LogOrigin("GJonesSpline", "getSplineVal()", WHERE));
1130 :
1131 : // Use GILDAS library routine, getbspl
1132 0 : Int numOfknots=knots.nelements();
1133 : Int failflag;
1134 : Bool dum;
1135 : Double yval;
1136 0 : getbspl(&numOfknots,
1137 : knots.getStorage(dum),
1138 : coeff.getStorage(dum),
1139 : &x,
1140 : &yval,
1141 : &failflag);
1142 0 : return yval;
1143 0 : }
1144 :
1145 : //----------------------------------------------------------------------------
1146 0 : void GJonesSpline::plotsolve(const Vector<Double>& x,
1147 : const Matrix<Double>& yall,
1148 : const Matrix<Double>& weightall,
1149 : const Vector<Double>& errall,
1150 : Matrix<Double>& coeff, Bool phasesoln){
1151 :
1152 :
1153 :
1154 : // Function to plot and compare Bandpass data and solution
1155 : // Input:
1156 : // x - vector of frequecies
1157 : // yall - matrix of data - shape is yall(freqindex, baselineindex)
1158 : // weightall - matrix of weight; same shape as yall
1159 : // errall - error of fits for each baseline
1160 : // coeff - matrix coefficients of polynomial fits -shape coeff(nant, degree)
1161 : // phasesoln - either phase or amplitude plot
1162 :
1163 :
1164 0 : LogIO os (LogOrigin("GJonesSpline", "plotsolve()", WHERE));
1165 :
1166 0 : String device;
1167 0 : device = calTableName();
1168 0 : if (phasesoln){
1169 0 : device = device + ".PHASE.ps/cps";
1170 : }
1171 : else{
1172 0 : device = device + ".AMP.ps/cps";
1173 : }
1174 :
1175 : os << LogIO::NORMAL
1176 : << "Generating plot file: " << device
1177 0 : << LogIO::POST;
1178 :
1179 0 : PGPlotter pg(device);
1180 0 : pg.subp(4,3);
1181 0 : Vector<Double> knots;
1182 0 : getKnots(x,knots);
1183 0 : Int numpoints= max(x.shape().asVector());
1184 0 : Int numplotpoints=20*numpoints;
1185 0 : Int num_valid_points=0;
1186 0 : Vector<Float> x1(numpoints);
1187 0 : Vector<Float> y1(numpoints);
1188 0 : Vector<Float> errarray(numpoints);
1189 0 : Vector<Double> xval(numplotpoints);
1190 0 : Vector<Float> xfloatval(numplotpoints);
1191 0 : Vector<Float> soly1(numplotpoints);
1192 0 : Vector<Double> y, weight;
1193 : Double err;
1194 0 : Int num_ant = nAnt();
1195 0 : Vector<Double> ant1coeff, ant2coeff;
1196 :
1197 : // Float max_data, min_data, max_err;
1198 : // max_err=max(errall);
1199 : // max_data=max(yall)+1.5*max_err;
1200 : // min_data=min(yall)-1.5*max_err;
1201 :
1202 0 : for (Int ant1=0; ant1 < num_ant; ++ant1){
1203 0 : for (Int ant2=ant1+1; ant2 < num_ant; ++ant2){
1204 0 : Int basnum=ant1*num_ant-ant1*(ant1+1)/2+ant2-1-ant1;
1205 :
1206 0 : x1.resize(numpoints);
1207 0 : y1.resize(numpoints);
1208 0 : errarray.resize(numpoints);
1209 0 : y.resize();
1210 0 : weight.resize();
1211 0 : ant1coeff=coeff.row(ant1);
1212 0 : ant2coeff=coeff.row(ant2);
1213 0 : y=yall.column(basnum);
1214 0 : weight=weightall.column(basnum);
1215 0 : err=errall[basnum];
1216 0 : num_valid_points=0;
1217 :
1218 0 : for(Int k=0; k < numpoints; k++){
1219 0 : if (weight(k)>0){
1220 0 : x1(num_valid_points)=x(k)-x(0);
1221 0 : y1(num_valid_points)=y(k);
1222 0 : if(phasesoln==true){
1223 0 : y1(num_valid_points)=remainder(y1(num_valid_points), 2*C::pi);
1224 0 : y1(num_valid_points)=y1(num_valid_points)/C::pi*180.0;
1225 0 : errarray(num_valid_points)=Float(err)*180.0/C::pi;
1226 : }
1227 : else{
1228 0 : y1(num_valid_points)=exp(y1(num_valid_points));
1229 0 : errarray(num_valid_points)=Float(err)*y1(num_valid_points);
1230 : }
1231 0 : num_valid_points+=1;
1232 : }; // weight(k) > 0
1233 : }; // for k=0
1234 :
1235 0 : if(num_valid_points == 0){
1236 0 : cout << "No valid point on baselines with antennas "
1237 0 : << ant1 << " & " << ant2 << endl;
1238 0 : return;
1239 : }
1240 :
1241 0 : for(Int k=0; k < numplotpoints; k++){
1242 0 : xval[k]= (k)*((x[numpoints-1]-x[0]))/Double(numplotpoints-1)
1243 0 : +(x[0]);
1244 :
1245 0 : if(phasesoln==true){
1246 0 : soly1(k)=0;
1247 : Double xx, yy1, yy2;
1248 0 : xx=xval[k];
1249 0 : yy1=getSplineVal(xx, knots, ant1coeff);
1250 0 : yy2=getSplineVal(xx, knots, ant2coeff);
1251 0 : soly1(k)=yy2-yy1;
1252 0 : while(soly1(k) > C::pi){
1253 0 : soly1(k) -= 2.0*(C::pi);
1254 : }
1255 0 : while(soly1(k) < (-(C::pi))){
1256 0 : soly1(k) += 2.0*(C::pi);
1257 : }
1258 0 : soly1(k)=soly1(k)/(C::pi)*180.0;
1259 :
1260 : }
1261 : else{
1262 :
1263 : Double xx, yy1, yy2;
1264 0 : xx=xval[k];
1265 0 : yy1=getSplineVal(xx, knots, ant1coeff);
1266 0 : yy2=getSplineVal(xx, knots, ant2coeff);
1267 0 : soly1[k]=(yy1+yy2);
1268 0 : soly1[k]=exp(soly1[k]);
1269 :
1270 :
1271 :
1272 : }
1273 :
1274 0 : xfloatval(k)=Float(xval(k)-xval(0));
1275 : }
1276 :
1277 :
1278 : //cout << " Mean of sol " << mean(soly1) << endl;
1279 0 : x1.resize(num_valid_points, true);
1280 0 : y1.resize(num_valid_points, true);
1281 0 : errarray.resize(num_valid_points, true);
1282 0 : pg.sci(1);
1283 :
1284 : Float max_data, min_data, max_err;
1285 0 : max_err=max(errarray);
1286 0 : max_data=max(y1);
1287 0 : min_data=min(y1);
1288 0 : max_data=max(max_data, max(soly1))+1.5*max_err;
1289 0 : min_data=min(min_data, min(soly1))-1.5*max_err;
1290 :
1291 0 : Float min_x= min(xfloatval);
1292 0 : Float max_x= max(xfloatval);
1293 :
1294 0 : Float edge=(max_x -min_x)*0.05;
1295 0 : min_x-=edge;
1296 0 : max_x+=edge;
1297 :
1298 0 : pg.sch(1.5);
1299 : // pg.env(min_x, max_x, min_data, max_data, 0, 0);
1300 0 : pg.page();
1301 0 : pg.svp(0.1,0.9,0.15,0.9);
1302 0 : pg.swin(min_x,max_x,min_data,max_data);
1303 0 : pg.box("BCNST",0.0,0,"BCNST",0.0,0);
1304 :
1305 0 : ostringstream oos;
1306 0 : if(phasesoln){
1307 0 : oos << "G SPLINE phase for baseline " << ant1+1 << " & " << ant2+1 ;
1308 0 : pg.lab("Time in s", "Phase in deg", " ");
1309 0 : pg.mtxt("T",0.75,0.5,0.5,oos);
1310 : }
1311 : else{
1312 0 : oos << "G SPLINE amplitude for baseline " << ant1+1 << " & " << ant2+1 ;
1313 0 : pg.lab("Time in s", "Amplitude", " ");
1314 0 : pg.mtxt("T",0.75,0.5,0.5,oos);
1315 : }
1316 0 : pg.sch(0.75);
1317 0 : pg.sci(2);
1318 0 : pg.errb(6, x1, y1, errarray, 2.0);
1319 0 : pg.sci(3);
1320 0 : pg.pt(x1, y1, 17);
1321 0 : pg.sci(4);
1322 0 : pg.line(xfloatval, soly1);
1323 :
1324 0 : }
1325 : }
1326 :
1327 0 : }
1328 :
1329 : //---------------------------------------------------------------------------
1330 :
1331 0 : Vector<Int> GJonesSpline::fieldIdRange()
1332 : {
1333 : // Return all field id.'s in the underlying MS
1334 : // Output:
1335 : // fieldIdRange Vector<Int> All FIELD_ID's in the MS
1336 : //
1337 : // Open the FIELD sub-table
1338 :
1339 0 : if (!vs_p) throw(AipsError("Error in GJonesSpline::fieldIdRange()"));
1340 :
1341 0 : const MSColumns& mscol(vs_p->iter().msColumns());
1342 0 : const MSFieldColumns& fldCol(mscol.field());
1343 :
1344 : // Fill vector containing all field id.'s
1345 0 : Vector<Int> retval(fldCol.nrow());
1346 0 : indgen(retval);
1347 :
1348 0 : return retval;
1349 0 : };
1350 :
1351 : //----------------------------------------------------------------------------
1352 :
1353 : /* TBD... (no raw phase processing, for now)
1354 :
1355 : void GJonesSpline::setRawPhaseVisSet(VisSet& vs){
1356 :
1357 : Block<Int> columns(0);
1358 :
1359 : rawvs_p = new VisSet (vs, columns, DBL_MAX);
1360 : rawPhaseRemoval_p=true;
1361 : fillRawPhaseBuff();
1362 : }
1363 :
1364 : void GJonesSpline::fillRawPhaseBuff(){
1365 :
1366 : //Make a phase array of BaselineIndex and TimeIndex of raw phases
1367 : LogIO os (LogOrigin("GJonesSpline", "fillRawPhaseBuff()", WHERE));
1368 :
1369 :
1370 : // Use the iterator from the underlying visibility set,
1371 : // and attach a visibility data buffer
1372 : VisIter& vi(rawvs_p->iter());
1373 : VisBuffer vb(vi);
1374 :
1375 : // Initialize stuff
1376 : // SimpleOrderedMap<String,Int> timeValueMap(0);
1377 : Vector<Double> timeValues;
1378 : PtrBlock<Vector<Complex>* > rawVisTimeSeries;
1379 : PtrBlock<Vector<Double>* > weightTimeSeries;
1380 : Int nTimeSeries = 0;
1381 :
1382 : // Initialize the baseline index
1383 : Int nAnt = vs_p->numberAnt();
1384 : Int nBasl = nAnt * (nAnt - 1) / 2;
1385 :
1386 :
1387 : // Iterate, accumulating the averaged visibility time-series
1388 : Int chunk;
1389 :
1390 : for (chunk=0, vi.originChunks(); vi.moreChunks(); vi.nextChunk(), chunk++) {
1391 : // Extract the current visibility buffer spectral window id.
1392 : // and number for frequency channels
1393 : Int spwid = vi.spectralWindow();
1394 : Int nChan = rawvs_p->numberChan()(spwid);
1395 : String fieldName = vi.fieldName();
1396 :
1397 : os << LogIO::NORMAL << "Slot= " << (chunk+1) << ", field= "
1398 : << fieldName << ", spw= " << spwid+1
1399 : << ", nchan= " << nChan << LogIO::POST;
1400 :
1401 : //average the buffer across frequency
1402 : for(vi.origin(); vi.more(); vi++){
1403 :
1404 : vb.visibility()=vb.correctedVisibility();
1405 : vb.freqAverage();
1406 :
1407 : // cout << "AFTER freq aver length" << vb.correctedVisibility().shape().asVector() << endl;
1408 :
1409 : // Iterate of the current visibility buffer, accumulating
1410 : // a common time series in the visibility data
1411 : for (Int row=0; row < vb.nRow(); row++) {
1412 : // Antenna numbers
1413 : Int ant1num = vb.antenna1()(row);
1414 : Int ant2num = vb.antenna2()(row);
1415 :
1416 : // Reject auto-correlation data
1417 : if (vb.antenna1()(row) != vb.antenna2()(row)) {
1418 : // Compute baseline index
1419 : Int baselineIndex = ant1num * nAnt - ant1num * (ant1num + 1) / 2 +
1420 : ant2num - 1 - ant1num;
1421 :
1422 : // Weight
1423 : Double rowWeight = vb.weight()(row);
1424 : if (rowWeight > 0) {
1425 : // Map the current time stamp to a time series index
1426 : // using 0.1 second precision, i.e. time values within
1427 : // a tenth of a second of each other will be accumulated
1428 : // within the same bin.
1429 : MVTime mvt(vb.time()(row) / C::day);
1430 : String timeKey = mvt.string(MVTime::TIME, 7);
1431 : // cout << "FILL Bef timekey "<< timeKey << " ntimeseries " << nTimeSeries << endl;
1432 : Int timeIndex = 0;
1433 : // Check the time stamp index to this precision
1434 : if (timeValueMap_p.isDefined(timeKey)) {
1435 : timeIndex = timeValueMap_p(timeKey);
1436 : } else {
1437 : // Create a new time series entry
1438 : timeIndex = nTimeSeries++;
1439 : // cout << "FILL timekey "<< timeKey << " index " << timeIndex << endl;
1440 : timeValueMap_p.define(timeKey, timeIndex);
1441 : timeValues.resize(nTimeSeries, true);
1442 : timeValues(timeIndex) = vb.time()(row);
1443 : rawVisTimeSeries.resize(nTimeSeries, true);
1444 : Complex czero(0,0);
1445 : rawVisTimeSeries[timeIndex] = new Vector<Complex>(nBasl, czero);
1446 : weightTimeSeries.resize(nTimeSeries, true);
1447 : weightTimeSeries[timeIndex] = new Vector<Double>(nBasl, 0);
1448 : };
1449 :
1450 : //CAUTION going to use corrected
1451 : // Accumulate the current visbility row in the common time series
1452 : (*rawVisTimeSeries[timeIndex])(baselineIndex) +=
1453 : vb.visibility()(0,row)(0) * rowWeight;
1454 : (*weightTimeSeries[timeIndex])(baselineIndex) += rowWeight;
1455 : };
1456 : };
1457 : };
1458 : };
1459 : }; // for (chunk...) iteration
1460 :
1461 : Int nTimes = timeValueMap_p.ndefined();
1462 : rawPhase_p.resize(nTimes, nBasl);
1463 : for (Int k=0; k< nTimes; ++k){
1464 : for(Int j=0; j< nBasl; ++j){
1465 : //rawPhase_p(k, j)=arg((*(rawVisTimeSeries[k]))(j));
1466 : rawPhase_p(k,j)=arg(((*rawVisTimeSeries[k])(j))
1467 : / ((*weightTimeSeries[k])(j)));
1468 :
1469 : }
1470 :
1471 : }
1472 :
1473 :
1474 :
1475 : }
1476 :
1477 :
1478 : Double GJonesSpline::getRawPhase(Int ant1, Int ant2, Double time){
1479 :
1480 : if (ant1 == ant2) return 0.0;
1481 : Int nAnt = vs_p->numberAnt();
1482 : Bool flip=false;
1483 : if(ant1 > ant2){
1484 : Int tmp=ant1;
1485 : ant1=ant2;
1486 : ant2=tmp;
1487 : flip=true;
1488 : }
1489 : Int baselineIndex = ant1 * nAnt - ant1 * (ant1 + 1) / 2 +
1490 : ant2 - 1 - ant1;
1491 :
1492 : MVTime mvt(time / C::day);
1493 : String timeKey = mvt.string(MVTime::TIME, 7);
1494 : Int timeIndex=timeValueMap_p(timeKey);
1495 : // if (ant1==0 && ant2==1)
1496 : // cout << "phase 1 2 " << rawPhase_p(timeIndex, baselineIndex)*180.0/C::pi << " Time " << timeKey << endl;
1497 : // cout << " Timekey " << timeKey << "Timeindex " << timeIndex << endl;
1498 : if(flip){
1499 : return -rawPhase_p(timeIndex, baselineIndex);
1500 : }
1501 : else {
1502 : return rawPhase_p(timeIndex, baselineIndex);
1503 : }
1504 :
1505 : }
1506 :
1507 : */
1508 :
1509 0 : void GJonesSpline::writeAsciiLog(const String& filename, const Matrix<Double>& coeff, const Vector<Double>& rmsFit, Bool phasesoln){
1510 :
1511 0 : Int numAnt=coeff.shape().asVector()(0);
1512 0 : const char* fil=filename.chars();
1513 0 : ofstream outLog(fil);
1514 :
1515 0 : if(phasesoln){
1516 0 : outLog << "*****PHASE summary********"<< endl;
1517 : }
1518 : else{
1519 0 : outLog << "*****AMPLITUDE summary********"<< endl;
1520 : }
1521 0 : outLog << "Antenna based spline coefficients " << endl;
1522 0 : for (Int k=0; k < numAnt; k++){
1523 0 : outLog << " " << k+1 << " " << coeff.row(k) << endl;
1524 :
1525 : }
1526 :
1527 :
1528 0 : outLog << " RMS Fit per baseline " << endl;
1529 0 : outLog << "Antenna1 Antenna2 rmsfit" << endl;
1530 0 : for (Int ant1=0; ant1 < numAnt; ++ant1){
1531 0 : for (Int ant2=ant1+1; ant2 < numAnt; ++ant2){
1532 0 : Int basnum=ant1*numAnt-ant1*(ant1+1)/2+ant2-1-ant1;
1533 0 : outLog << " " << ant1+1 << " " << ant2+1 << " "
1534 0 : << rmsFit[basnum] << endl;
1535 : }
1536 : }
1537 0 : }
1538 :
1539 : } //# NAMESPACE CASA - END
1540 :
|
