Line data    Source code

       1             : //# AccorJones.cc: Implementation of AccorJones
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011,2017
       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: aips2-request@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/AccorJones.h>
      28             : 
      29             : #include <msvis/MSVis/VisBuffer.h>
      30             : #include <msvis/MSVis/VisBuffAccumulator.h>
      31             : #include <synthesis/MeasurementEquations/VisEquation.h>
      32             : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
      33             : 
      34             : #include <casacore/casa/Arrays/ArrayMath.h>
      35             : #include <casacore/casa/Arrays/MatrixMath.h>
      36             : #include <casacore/casa/BasicSL/String.h>
      37             : #include <casacore/casa/Utilities/Assert.h>
      38             : #include <casacore/casa/Exceptions/Error.h>
      39             : #include <casacore/casa/System/Aipsrc.h>
      40             : 
      41             : #include <sstream>
      42             : 
      43             : 
      44             : 
      45             : using namespace casacore;
      46             : namespace casa { //# NAMESPACE CASA - BEGIN
      47             : 
      48             : // **********************************************************
      49             : //  AccorJones Implementations
      50             : //
      51             : 
      52           0 : AccorJones::AccorJones(VisSet& vs) :
      53             :   VisCal(vs),             // virtual base
      54             :   VisMueller(vs),         // virtual base
      55           0 :   SolvableVisJones(vs)    // immediate parent
      56             : {
      57           0 :   if (prtlev()>2) cout << "Accor::Accor(vs)" << endl;
      58           0 : }
      59             : 
      60           0 : AccorJones::AccorJones(String msname,Int MSnAnt,Int MSnSpw) :
      61             :   VisCal(msname,MSnAnt,MSnSpw),             // virtual base
      62             :   VisMueller(msname,MSnAnt,MSnSpw),         // virtual base
      63           0 :   SolvableVisJones(msname,MSnAnt,MSnSpw)    // immediate parent
      64             : {
      65           0 :   if (prtlev()>2) cout << "Accor::Accor(msname,MSnAnt,MSnSpw)" << endl;
      66           0 : }
      67             : 
      68           0 : AccorJones::AccorJones(const MSMetaInfoForCal& msmc) :
      69             :   VisCal(msmc),             // virtual base
      70             :   VisMueller(msmc),         // virtual base
      71           0 :   SolvableVisJones(msmc)    // immediate parent
      72             : {
      73           0 :   if (prtlev()>2) cout << "Accor::Accor(msmc)" << endl;
      74           0 : }
      75             : 
      76           0 : AccorJones::AccorJones(const Int& nAnt) :
      77             :   VisCal(nAnt),
      78             :   VisMueller(nAnt),
      79           0 :   SolvableVisJones(nAnt)
      80             : {
      81           0 :   if (prtlev()>2) cout << "Accor::Accor(nAnt)" << endl;
      82           0 : }
      83             : 
      84           0 : AccorJones::~AccorJones() {
      85           0 :   if (prtlev()>2) cout << "Accor::~Accor()" << endl;
      86           0 : }
      87             : 
      88           0 : void AccorJones::selfSolveOne(VisBuffGroupAcc&) {
      89             : 
      90             :   // NB: Use "AccorJones::selfSolveOne(SDBList& sdbs)" instead!
      91           0 :   throw(AipsError("VisBuffGroupAcc is invalid for AccorJones"));
      92             : }
      93             : 
      94           0 : void AccorJones::selfSolveOne(SDBList& sdbs) {
      95             : 
      96             :   // Use just the first SDB in the SDBList, for now
      97           0 :   SolveDataBuffer& sdb(sdbs(0));
      98             : 
      99             :   // Make a guess at antenna-based G
     100             :   //  Correlation membership-dependencexm
     101             :   //  assumes corrs in canonical order
     102             :   //  nCorr = 1: use icorr=0
     103             :   //  nCorr = 2: use icorr=[0,1]
     104             :   //  nCorr = 4: use icorr=[0,3]
     105             : 
     106           0 :   Int nCorr(2);
     107           0 :   Int nDataCorr(sdb.nCorrelations());
     108           0 :   Vector<Int> corridx(nCorr,0);
     109           0 :   if (nDataCorr==2) {
     110           0 :     corridx(0)=0;
     111           0 :     corridx(1)=1;
     112             :   }
     113           0 :   else if (nDataCorr==4) {
     114           0 :     corridx(0)=0;
     115           0 :     corridx(1)=3;
     116             :   }
     117             : 
     118           0 :   Int guesschan(sdb.nChannels()-1);
     119             : 
     120             :   //  cout << "guesschan = " << guesschan << endl;
     121             :   //  cout << "nCorr = " << nCorr << endl;
     122             :   //  cout << "corridx = " << corridx << endl;
     123             : 
     124             :   // Find out which ants are available
     125           0 :   Int nRow=sdb.nRows();
     126           0 :   Vector<Bool> rowok(nRow,False);
     127           0 :   for (Int irow=0;irow<nRow;++irow) {
     128             : 
     129           0 :     const Int& a1(sdb.antenna1()(irow));
     130           0 :     const Int& a2(sdb.antenna2()(irow));
     131             : 
     132             :     // Is this row ok?
     133           0 :     rowok(irow)= (!sdb.flagRow()(irow) &&
     134           0 :                   a1==a2);
     135             :   }
     136             : 
     137           0 :   const Cube<Complex>& V(sdb.visCubeCorrected());
     138           0 :   Float amp(0.0);
     139           0 :   solveCPar()=Complex(1.0);
     140           0 :   solveParOK()=False;
     141           0 :   for (Int irow=0;irow<nRow;++irow) {
     142             : 
     143           0 :     if (rowok(irow)) {
     144           0 :       const Int& a1=sdb.antenna1()(irow);
     145           0 :       const Int& a2=sdb.antenna2()(irow);
     146             : 
     147           0 :       if (a1 == a2) {
     148             : 
     149           0 :         for (Int icorr=0;icorr<nCorr;icorr++) {
     150           0 :           const Complex& Vi(V(corridx(icorr),guesschan,irow));
     151           0 :           amp=abs(Vi);
     152           0 :           if (amp>0.0f) {
     153           0 :             solveCPar()(icorr,0,a1) = sqrt(amp);
     154           0 :             solveParOK()(icorr,0,a1) = True;
     155             :           }
     156             :         }
     157             :       }
     158             :     } // rowok
     159             :   } // irow
     160             : 
     161             :   // For scalar data, set "other" pol soln to zero
     162           0 :   if (nDataCorr == 1)
     163           0 :     solveCPar()(IPosition(3,1,0,0),IPosition(3,1,0,nAnt()-1))=Complex(0.0);
     164           0 : }
     165             : 
     166             : // Fill the trivial DJ matrix elements
     167           0 : void AccorJones::initTrivDJ() {
     168             : 
     169           0 :   if (prtlev()>4) cout << "   Accor::initTrivDJ" << endl;
     170             : 
     171             :   // Must be trivial
     172           0 :   AlwaysAssert((trivialDJ()),AipsError);
     173             : 
     174             :   //  1 0     0 0
     175             :   //  0 0  &  0 1
     176             :   //
     177           0 :   if (diffJElem().nelements()==0) {
     178           0 :     diffJElem().resize(IPosition(4,2,2,1,1));
     179           0 :     diffJElem()=0.0;
     180           0 :     diffJElem()(IPosition(4,0,0,0,0))=Complex(1.0);
     181           0 :     diffJElem()(IPosition(4,1,1,0,0))=Complex(1.0);
     182             :   }
     183             : 
     184           0 : }
     185             : 
     186           0 : void AccorJones::keepNCT() {
     187           0 :   for (Int elem=0;elem<nElem();++elem) {
     188           0 :     ct_->addRow(1);
     189             : 
     190           0 :     CTMainColumns ncmc(*ct_);
     191             : 
     192             :     // We are adding to the most-recently added row
     193           0 :     Int row=ct_->nrow()-1;
     194             : 
     195             :     // Meta-info
     196           0 :     ncmc.time().put(row,refTime());
     197           0 :     ncmc.fieldId().put(row,currField());
     198           0 :     ncmc.spwId().put(row,currSpw());
     199           0 :     ncmc.scanNo().put(row,currScan());
     200           0 :     ncmc.obsId().put(row,currObs());
     201           0 :     ncmc.interval().put(row,0.0);
     202           0 :     ncmc.antenna1().put(row,elem);
     203           0 :     ncmc.antenna2().put(row,-1);
     204             : 
     205             :     // Params
     206           0 :     ncmc.cparam().put(row,solveAllCPar().xyPlane(elem));
     207             : 
     208             :     // Stats
     209           0 :     ncmc.paramerr().put(row,solveAllParErr().xyPlane(elem));
     210           0 :     ncmc.snr().put(row,solveAllParErr().xyPlane(elem));
     211           0 :     ncmc.flag().put(row,!solveAllParOK().xyPlane(elem));
     212             :   }
     213             : 
     214             :   // This spw now has some solutions in it
     215           0 :   spwOK_(currSpw())=True;
     216           0 : }
     217             : 
     218             : } //# NAMESPACE CASA - END