diff --git a/CEP/DP3/DPPP/include/DPPP/phasefitter.h b/CEP/DP3/DPPP/include/DPPP/phasefitter.h index 82bff267fc166ab8f4ed2a8887835baa581625ac..f09a2db8eb6ee017ee32b71a2a34f458ea58921c 100644 --- a/CEP/DP3/DPPP/include/DPPP/phasefitter.h +++ b/CEP/DP3/DPPP/include/DPPP/phasefitter.h @@ -285,6 +285,7 @@ class PhaseFitter void bruteForceSearchTEC2Model(double& lowerAlpha, double& upperAlpha, double& beta) const; double ternarySearchTEC2ModelAlpha(double startAlpha, double endAlpha, double& beta) const; void fillDataWithTEC2Model(double alpha, double beta); + void fillDataWithTEC1Model(double alpha); void bruteForceSearchTEC1Model(double& lowerAlpha, double& upperAlpha) const; double ternarySearchTEC1ModelAlpha(double startAlpha, double endAlpha) const; diff --git a/CEP/DP3/DPPP/src/GainCal.cc b/CEP/DP3/DPPP/src/GainCal.cc index e863212c633ba7122471df3d18c3f5715a076b7a..8bfb2148014cbed5ac984842ffceb421ec51c1df 100644 --- a/CEP/DP3/DPPP/src/GainCal.cc +++ b/CEP/DP3/DPPP/src/GainCal.cc @@ -124,7 +124,7 @@ namespace LOFAR { ASSERT(itsMode=="diagonal" || itsMode=="phaseonly" || itsMode=="fulljones" || itsMode=="scalarphase" || itsMode=="amplitudeonly" || itsMode=="scalaramplitude" || - itsMode=="tec"); + itsMode=="tecandphase" || itsMode=="tec"); } GainCal::~GainCal() @@ -172,7 +172,7 @@ namespace LOFAR { itsSols.reserve(itsTimeSlotsPerParmUpdate); // Initialize phase fitters, set their frequency data - if (itsMode=="tec") { + if (itsMode=="tecandphase" || itsMode=="tec") { itsTECSols.reserve(itsTimeSlotsPerParmUpdate); itsPhaseFitters.reserve(itsNFreqCells); // TODO: could be numthreads instead @@ -260,7 +260,7 @@ namespace LOFAR { FlagCounter::showPerc1 (os, itsTimerSolve.getElapsed(), totaltime); os << " of it spent in estimating gains and computing residuals" << endl; - if (itsMode == "tec") { + if (itsMode == "tec" || itsMode == "tecandphase") { os << " "; FlagCounter::showPerc1 (os, itsTimerPhaseFit.getElapsed(), totaltime); os << " of it spent in fitting phases" << endl; @@ -451,7 +451,7 @@ namespace LOFAR { continue; } - if (itsMode=="tec") { + if (itsMode=="tec" || itsMode=="tecandphase") { iS[ch/itsNChan].incrementWeight(weight[bl*nCr*nCh+ch*nCr]); } @@ -540,7 +540,7 @@ namespace LOFAR { uint iter=0; - casa::Matrix<double> tecsol(2, info().antennaNames().size(), 0); + casa::Matrix<double> tecsol(itsMode=="tecandphase"?2:1, info().antennaNames().size(), 0); std::vector<StefCal::Status> converged(itsNFreqCells,StefCal::NOTCONVERGED); for (;iter<itsMaxIter;++iter) { @@ -556,7 +556,7 @@ namespace LOFAR { } } - if (itsMode=="tec") { + if (itsMode=="tec" || itsMode=="tecandphase") { itsTimerSolve.stop(); itsTimerPhaseFit.start(); casa::Matrix<casa::DComplex> sols_f(itsNFreqCells, info().antennaNames().size()); @@ -624,7 +624,11 @@ namespace LOFAR { if (numpoints>1) { // TODO: limit should be higher //cout<<"tecsol(0,"<<st<<")="<<tecsol(0,st)<<", tecsol(1,"<<st<<")="<<tecsol(1,st)<<endl; - cost=itsPhaseFitters[st]->FitDataToTEC2Model(tecsol(0, st), tecsol(1,st)); + if (itsMode=="tecandphase") { + cost=itsPhaseFitters[st]->FitDataToTEC2Model(tecsol(0, st), tecsol(1,st)); + } else { // itsMode=="tec" + cost=itsPhaseFitters[st]->FitDataToTEC1Model(tecsol(0, st)); + } // Update solution in stefcal for (uint freqCell=0; freqCell<itsNFreqCells; ++freqCell) { ASSERT(isFinite(phases[freqCell])); @@ -632,7 +636,9 @@ namespace LOFAR { } } else { tecsol(0, st) = 0; //std::numeric_limits<double>::quiet_NaN(); - tecsol(1, st) = 0; //std::numeric_limits<double>::quiet_NaN(); + if (itsMode=="tecandphase") { + tecsol(1, st) = 0; //std::numeric_limits<double>::quiet_NaN(); + } } if (st==34 && itsDebugLevel>0) { @@ -706,7 +712,7 @@ namespace LOFAR { } } itsSols.push_back(sol); - if (itsMode=="tec") { + if (itsMode=="tec" || itsMode=="tecandphase") { itsTECSols.push_back(tecsol); } @@ -793,7 +799,7 @@ namespace LOFAR { name=string("CommonScalarPhase:"); } else if (itsMode=="scalaramplitude") { name=string("CommonScalarAmplitude:"); - } else if (itsMode=="tec") { + } else if (itsMode=="tec" || itsMode=="tecandphase") { name=string("TEC:"); } else { @@ -881,9 +887,9 @@ namespace LOFAR { itsMode=="tec") && pol>0) { continue; } - int realimmax; // For tec, this functions as dummy between tec and commonscalarphase + int realimmax; // For tecandphase, this functions as dummy between tec and commonscalarphase if (itsMode=="phaseonly" || itsMode=="scalarphase" || - itsMode=="amplitudeonly" || itsMode=="scalaramplitude") { + itsMode=="amplitudeonly" || itsMode=="scalaramplitude" || itsMode=="tec") { realimmax=1; } else { realimmax=2; @@ -901,7 +907,7 @@ namespace LOFAR { name=name+strri[realim]; } } - if (itsMode=="tec") { + if (itsMode=="tecandphase" || itsMode=="tec") { if (realim==0) { name="TEC:"; } else { @@ -930,7 +936,7 @@ namespace LOFAR { values(freqCell, ts) = arg(itsSols[ts](pol/3,st,freqCell)); } else if (itsMode=="scalaramplitude" || itsMode=="amplitudeonly") { values(freqCell, ts) = abs(itsSols[ts](pol/3,st,freqCell)); - } else if (itsMode=="tec") { + } else if (itsMode=="tec" || itsMode=="tecandphase") { if (realim==0) { values(freqCell, ts) = itsTECSols[ts](realim,st) / 8.44797245e9; } else { diff --git a/CEP/DP3/DPPP/src/phasefitter.cc b/CEP/DP3/DPPP/src/phasefitter.cc index 07c171dd743654cd1cbf7455598bbe0138650894..2391ada19965fd23ca752007fd09718170ade858 100644 --- a/CEP/DP3/DPPP/src/phasefitter.cc +++ b/CEP/DP3/DPPP/src/phasefitter.cc @@ -114,6 +114,12 @@ void PhaseFitter::fillDataWithTEC2Model(double alpha, double beta) _phases[ch] = TEC2ModelFunc(_frequencies[ch], alpha, beta); } +void PhaseFitter::fillDataWithTEC1Model(double alpha) +{ + for(size_t ch=0; ch!=Size(); ++ch) + _phases[ch] = TEC1ModelFuncWrapped(_frequencies[ch], alpha); +} + void PhaseFitter::FitTEC2ModelParameters(double& alpha, double& beta) const { double lowerAlpha = -40000.0e6, upperAlpha = 40000.0e6; @@ -130,6 +136,13 @@ double PhaseFitter::FitDataToTEC2Model(double& alpha, double& beta) return TEC2ModelCost(alpha, beta); } +double PhaseFitter::FitDataToTEC1Model(double& alpha) +{ + FitTEC1ModelParameters(alpha); + fillDataWithTEC1Model(alpha); + return TEC1ModelCost(alpha); +} + void PhaseFitter::FitTEC1ModelParameters(double& alpha) const { double lowerAlpha = -40000.0e6, upperAlpha = 40000.0e6;