LCOV - code coverage report
Current view: top level - synthesis/MeasurementComponents - AMueller.cc (source / functions) Hit Total Coverage
Test: ctest_coverage.info Lines: 211 274 77.0 %
Date: 2023-11-06 10:06:49 Functions: 19 36 52.8 %

          Line data    Source code
       1             : //# AMueller.cc: Implementation of AMueller
       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: 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/CalTables/CalDescColumns.h>
      28             : #include <synthesis/CalTables/CTColumns.h>
      29             : #include <synthesis/CalTables/CTMainColumns.h>
      30             : //#include <ms/MeasurementSets/MSColumns.h>
      31             : #include <msvis/MSVis/VisBuffer.h>
      32             : #include <msvis/MSVis/VisBuffGroupAcc.h>
      33             : #include <msvis/MSVis/VBContinuumSubtractor.h>
      34             : #include <synthesis/MeasurementComponents/CalCorruptor.h>
      35             : #include <synthesis/MeasurementComponents/AMueller.h>
      36             : #include <synthesis/MeasurementEquations/VisEquation.h>
      37             : 
      38             : #include <msvis/MSVis/VisBuffer2.h>
      39             : 
      40             : 
      41             : using namespace casacore;
      42             : namespace casa { //# NAMESPACE CASA - BEGIN
      43             : 
      44             : 
      45             : // **********************************************************
      46             : //  AMueller
      47             : //
      48             : 
      49             : 
      50          14 : AMueller::AMueller(VisSet& vs) :
      51             :   VisCal(vs),             // virtual base
      52             :   VisMueller(vs),         // virtual base
      53             :   MMueller(vs),            // immediate parent
      54             :   fitorder_p(0),
      55             :   doSub_p(true),
      56          14 :   nCorr_p(-1)
      57             : {
      58          14 :   if (prtlev()>2) cout << "A::A(vs)" << endl;
      59             : 
      60          14 :   init();
      61          14 : }
      62             : 
      63           0 : AMueller::AMueller(String msname,Int MSnAnt,Int MSnSpw) :
      64             :   VisCal(msname,MSnAnt,MSnSpw),             // virtual base
      65             :   VisMueller(msname,MSnAnt,MSnSpw),         // virtual base
      66             :   MMueller(msname,MSnAnt,MSnSpw),            // immediate parent
      67             :   fitorder_p(0),
      68             :   doSub_p(true),
      69           0 :   nCorr_p(-1)
      70             : {
      71           0 :   if (prtlev()>2) cout << "A::A(msname,MSnAnt,MSnSpw)" << endl;
      72             : 
      73           0 :   init();
      74           0 : }
      75             : 
      76             : 
      77           0 : AMueller::AMueller(const MSMetaInfoForCal& msmc) :
      78             :   VisCal(msmc),             // virtual base
      79             :   VisMueller(msmc),         // virtual base
      80             :   MMueller(msmc),            // immediate parent
      81             :   fitorder_p(0),
      82             :   doSub_p(true),
      83           0 :   nCorr_p(-1)
      84             : {
      85           0 :   if (prtlev()>2) cout << "A::A(msmc)" << endl;
      86             : 
      87           0 :   init();
      88           0 : }
      89             : 
      90             : 
      91           0 : AMueller::AMueller(const Int& nAnt) :
      92             :   VisCal(nAnt),
      93             :   VisMueller(nAnt),
      94             :   MMueller(nAnt),            // immediate parent
      95             :   fitorder_p(0),
      96             :   doSub_p(true),
      97           0 :   nCorr_p(-1)
      98             : {
      99           0 :   if (prtlev()>2) cout << "A::A(nAnt)" << endl;
     100             : 
     101           0 :   init();
     102           0 : }
     103             : 
     104          28 : AMueller::~AMueller() {
     105          14 :   if (prtlev()>2) cout << "A::~A()" << endl;
     106          28 : }
     107             : 
     108          14 : void AMueller::init()
     109             : {
     110          14 :   lofreq_p.resize(nSpw());
     111          14 :   hifreq_p.resize(nSpw());
     112          14 :   totnumchan_p.resize(nSpw());
     113          14 :   spwApplied_p.resize(nSpw());
     114             :   
     115          14 :   lofreq_p = -1.0;
     116          14 :   hifreq_p = -1.0;
     117          14 :   totnumchan_p = 0;
     118          14 :   spwApplied_p = false;
     119          14 : }
     120             : 
     121           7 : void AMueller::setSolve(const Record& solvepar) {
     122             : 
     123             :   // Call parent
     124           7 :   MMueller::setSolve(solvepar);
     125             : 
     126             :   // Extract the AMueller-specific solve parameters
     127           7 :   if(solvepar.isDefined("fitorder"))
     128           7 :     fitorder_p = solvepar.asInt("fitorder");
     129             : 
     130           7 :   nChanParList() = fitorder_p + 1;  // Orders masquerade as output chans.
     131             : 
     132             :   // Override preavg 
     133             :   // (solver will fail if we don't average completely in each solint)
     134           7 :   preavg()=DBL_MAX;
     135           7 : }
     136             : 
     137           7 : void AMueller::setSolveChannelization(VisSet& vs)
     138             : {
     139          14 :   Vector<Int> startDatChan(vs.startChan());
     140             : 
     141             :   // If fitorder_p != 0, this is frequency dependent.
     142           7 :   if(fitorder_p){
     143             :     // AMueller keeps its polynomial orders where channels would normally go,
     144             :     // and typically (advisedly) the number of orders is << the number of data
     145             :     // channels.  *Otherwise* the overall par shape follows the data shape.
     146           2 :     nChanParList() = fitorder_p + 1;  // Deja vu from setSolve().
     147           2 :     startChanList() = startDatChan;
     148             : 
     149             :     // However, for solving, we will only consider one channel at a time:
     150           2 :     nChanMatList() = 1;
     151             :   }
     152             :   else {
     153             :     // Pars are not themselves channel-dependent
     154           5 :     nChanParList() = 1;
     155             : 
     156             :     // Check if matrices may still be freq-dep:
     157           5 :     if (freqDepMat()) {
     158             :       // cal is an explicit f(freq) (e.g., like delay)
     159           0 :       nChanMatList()  = vs.numberChan();
     160           0 :       startChanList() = startDatChan;
     161             :     } else {
     162             :       // cal has no freq dep at all
     163           5 :       nChanMatList()  = Vector<Int>(nSpw(),1);
     164           5 :       startChanList() = Vector<Int>(nSpw(),0);
     165             :     }
     166             :   }
     167             : 
     168             :   // At this point:
     169             :   //  1. nChanParList() represents the number of coefficients per polynomial,
     170             :   //     appropriate for shaping the CalSet.
     171             :   //  2. nChanMatList() represents the per-Spw matrix channel axis length to
     172             :   //     be used during the solve, independent of the parameter channel
     173             :   //     axis length.  In the solve context, nChanMat()>1 when there is
     174             :   //     more than one channel of data upon which the (single channel)
     175             :   //     solve parameters depend (e.g. polynomial order != 1)
     176           7 : }
     177             : 
     178           7 : Int AMueller::sizeUpSolve(VisSet& vs, Vector<Int>& nChunkPerSol)
     179             : {
     180           7 :   sortVisSet(vs);
     181           7 :   VisIter& vi(vs.iter());
     182             : 
     183           7 :   nCorr_p = fitorder_p ? vi.nCorr() : 2;
     184             : 
     185             :   // Would SolvableVisCal be better here?  After all the need for this
     186             :   // specialization started with MMueller::nPar()...
     187           7 :   return MMueller::sizeUpSolve(vs, nChunkPerSol);
     188             : }
     189             : 
     190          39 : void AMueller::selfSolveOne(VisBuffGroupAcc& vbga)
     191             : {
     192             :   // Solver for the polynomial continuum fit.  It is overkill for fitorder_p ==
     193             :   // 0, and not used in that case.
     194             : 
     195         117 :   LogIO os(LogOrigin("AMueller", "selfSolveOne()", WHERE));
     196          78 :   VBContinuumSubtractor vbcs;
     197             : 
     198             :   // Initialize
     199          39 :   if(lofreq_p[currSpw()] < 0.0){  // 1st time for this spw, so let vbga
     200           2 :     vbcs.initFromVBGA(vbga);      // provide the info.
     201           2 :     lofreq_p[currSpw()] = vbcs.getLowFreq();
     202           2 :     hifreq_p[currSpw()] = vbcs.getHighFreq();
     203           2 :     totnumchan_p[currSpw()] = vbcs.getTotNumChan();
     204             :   }
     205             :   else                            // Reuse the prev vals for consistency.
     206          37 :     vbcs.init(solveParOK().shape(), nAnt() - 1, totnumchan_p[currSpw()],
     207          37 :               lofreq_p[currSpw()], hifreq_p[currSpw()]);
     208             : 
     209          39 :   vbcs.fit(vbga, fitorder_p, MS::DATA, solveCPar(), solveParOK(),
     210          39 :            false, false, !append());
     211          39 : }
     212             : 
     213           7 : void AMueller::storeNCT() {
     214             : 
     215             :   // Update the freq info
     216           7 :   if(fitorder_p != 0){    // Store lofreq_p[currSpw()] and hifreq_p[currSpw()]
     217             : 
     218             :     // Access SPW subtable
     219           4 :     CTColumns ncc(*ct_);
     220             : 
     221             :     // We are expecting nSpw rows...
     222           2 :     AlwaysAssert( Int(ncc.spectralWindow().chanFreq().nrow())==nSpw(), AipsError);
     223             : 
     224             :     // We store the freq domain as a pair of values in chanFreq(),
     225             :     //  and reset numChan()=fitorder_p
     226             : 
     227           4 :     Vector<Double> freqrange(2);
     228          25 :     for (Int ispw=0;ispw<nSpw();++ispw) {
     229             :       // Only if values are ok...
     230          25 :       if (lofreq_p[ispw]>0.0 &&
     231           2 :           hifreq_p[ispw]>0.0) {
     232           2 :         freqrange(0)=lofreq_p[ispw];
     233           2 :         freqrange(1)=hifreq_p[ispw];
     234           2 :         ncc.spectralWindow().chanFreq().setShape(ispw,IPosition(1,2));
     235           2 :         ncc.spectralWindow().chanFreq().put(ispw,freqrange);
     236           2 :         ncc.spectralWindow().numChan().put(ispw,fitorder_p+1);
     237           2 :         ncc.spectralWindow().flagRow().put(ispw,false);
     238             :       }
     239             :       else
     240             :         // Mark this spw flagged (no solutions in main
     241          21 :         ncc.spectralWindow().flagRow().put(ispw,true);
     242             :     }
     243             :   }
     244             : 
     245             :   // Now call conventional store
     246           7 :   MMueller::storeNCT();
     247             : 
     248           7 : }
     249             : 
     250           0 : void AMueller::hurl(const String& origin, const String& msg)
     251             : {
     252           0 :   LogIO os(LogOrigin("AMueller", origin, WHERE));
     253             :   
     254           0 :   os << msg << LogIO::EXCEPTION;
     255           0 : }
     256             : 
     257           7 : void AMueller::setApply(const Record& applypar)
     258             : {
     259          21 :   LogIO os(LogOrigin("AMueller", "setApply()", WHERE));
     260             :   
     261           7 :   if(applypar.isDefined("table")){
     262           7 :     calTableName() = applypar.asString("table");
     263           7 :     verifyCalTable(calTableName());
     264             : 
     265             :     {
     266             :       // TBD: Improve this kluge!
     267             :       //  e.g., fitorder from a keyword, etc.
     268          14 :       NewCalTable ct0(calTableName(),Table::Old,Table::Plain);
     269          14 :       ROCTMainColumns ncmc(ct0);
     270           7 :       IPosition sh=ncmc.cparam().shape(0);
     271           7 :       nCorr_p = sh(0);
     272           7 :       fitorder_p= sh(1)-1;
     273           7 :       nChanParList().set(sh(1));
     274             :     }
     275             :   }
     276             :   else
     277             :     os << "AMueller::setApply(Record&) needs the record to have a table"
     278           0 :        << LogIO::EXCEPTION;
     279             : 
     280             : 
     281             :   // Call parent to get the general stuff
     282           7 :   MMueller::setApply(applypar);
     283             : 
     284             :   // Now get some uvcontsub-specific stuff
     285             : 
     286             :   // Extract freq domain for non-trivial fits
     287           7 :   if(fitorder_p != 0){
     288             : 
     289             :     // Access to SPW subtable..
     290           4 :     ROCTColumns ncc(*ct_);
     291           4 :     Vector<Bool> spwflrow;
     292           2 :     ncc.spectralWindow().flagRow().getColumn(spwflrow);
     293           2 :     Int nspw=spwflrow.nelements();
     294             : 
     295             :     // TBD: What if caltable contains more or less than nSpw() (from MS) spws?
     296             : 
     297           2 :     lofreq_p.resize(nspw);
     298           2 :     hifreq_p.resize(nspw);
     299             : 
     300           4 :     Vector<Double> freqrange(2);
     301           2 :     if (spwMap()(0)==-999) {  // auto-fan-out
     302           1 :       Int ispw=0;
     303           1 :       while (spwflrow(ispw)) ++ispw;  // find first good spw
     304           1 :       ncc.spectralWindow().chanFreq().get(ispw,freqrange);
     305           1 :       lofreq_p.set(freqrange(0));
     306           1 :       hifreq_p.set(freqrange(1));
     307             :     }
     308             :     else {  // fill lo/hifreq per spw (only good ones)
     309           2 :       for (Int ispw=0;ispw<nspw;++ispw) {
     310           1 :         if (!spwflrow(ispw)) {
     311           1 :           ncc.spectralWindow().chanFreq().get(ispw,freqrange);
     312           1 :           lofreq_p[ispw]=freqrange(0);
     313           1 :           hifreq_p[ispw]=freqrange(1);
     314             :         }
     315             :       }
     316             :     }
     317             : 
     318             :   }
     319             : 
     320           7 : }
     321             : 
     322             : // Apply this calibration to VisBuffer visibilities
     323         102 : void AMueller::applyCal(VisBuffer& vb, Cube<Complex>& Vout,Bool trial)
     324             : {
     325         306 :   LogIO os(LogOrigin("AMueller", "applyCal()", WHERE));
     326             : 
     327         102 :   if(fitorder_p == 0){
     328         102 :     VisMueller::applyCal(vb, Vout, trial);
     329             :   }
     330             :   else{
     331           0 :     if(prtlev() > 3)
     332           0 :       os << "  AMueller::applyCal()" << LogIO::POST;
     333             : 
     334           0 :     if (trial)
     335           0 :       throw(AipsError("trial apply not supported by AMueller with fitorder_p>0"));
     336             :     
     337           0 :     Int cspw = currSpw();
     338           0 :     VBContinuumSubtractor vbcs;
     339             : 
     340           0 :     if (lofreq_p[cspw]>0.0 &&
     341           0 :         hifreq_p[cspw]>0.0)
     342           0 :       vbcs.init(currCPar().shape(), nAnt() - 1, totnumchan_p[cspw],
     343           0 :                 lofreq_p[cspw], hifreq_p[cspw]);
     344             :     else
     345           0 :       throw(AipsError("AMueller::applyCal: Bad freq domain info."));
     346             : 
     347             :     // Correct DATA (will write out to CORRECTED_DATA later)
     348           0 :     MS::PredefinedColumns whichcol = MS::DATA;
     349           0 :     if(!vbcs.apply(vb, whichcol, currCPar(), currParOK(), doSub_p,
     350           0 :                    !spwApplied_p[cspw]))
     351             :       throw(AipsError("Could not place the continuum-subtracted data in "
     352           0 :                       + MS::columnName(whichcol)));
     353           0 :     spwApplied_p[cspw] = true;
     354             :   }
     355         102 : }
     356             : 
     357             : 
     358             : // Apply this calibration to VisBuffer visibilities
     359         241 : void AMueller::applyCal2(vi::VisBuffer2& vb, 
     360             :                          casacore::Cube<casacore::Complex>& Vout,casacore::Cube<casacore::Float>& Wout,
     361             :                          casacore::Bool trial)
     362             : 
     363             : {
     364         723 :   LogIO os(LogOrigin("AMueller", "applyCal2()", WHERE));
     365             : 
     366         241 :   if(fitorder_p == 0){
     367         162 :     VisMueller::applyCal2(vb, Vout, Wout, trial);
     368             :   }
     369             :   else{
     370          79 :     if(prtlev() > 3)
     371           0 :       os << "  AMueller::applyCal2()" << LogIO::POST;
     372             : 
     373          79 :     if (trial)
     374           0 :       throw(AipsError("trial apply not supported by AMueller with fitorder_p>0"));
     375             :     
     376          79 :     Int cspw = currSpw();
     377          79 :     VBContinuumSubtractor vbcs;
     378             : 
     379         158 :     if (lofreq_p[cspw]>0.0 &&
     380          79 :         hifreq_p[cspw]>0.0)
     381          79 :       vbcs.init(currCPar().shape(), nAnt() - 1, totnumchan_p[cspw],
     382          79 :                 lofreq_p[cspw], hifreq_p[cspw]);
     383             :     else
     384           0 :       throw(AipsError("AMueller::applyCal2: Bad freq domain info."));
     385             : 
     386             :     // Correct Vout
     387          79 :     if(!vbcs.apply2(vb, Vout, currCPar(), currParOK(), doSub_p,
     388          79 :                    !spwApplied_p[cspw]))
     389           0 :       throw(AipsError("Error applying continuum-subtraction"));
     390             : 
     391          79 :     spwApplied_p[cspw] = true;
     392             :   }
     393         241 : }
     394             : 
     395         136 : void AMueller::corrupt(VisBuffer& vb)
     396             : {
     397         408 :   LogIO os(LogOrigin("AMueller", "corrupt()", WHERE));
     398             : 
     399         136 :   if(prtlev() > 3)
     400           0 :     os << LogIO::NORMAL << "  A::corrupt()" << LogIO::POST;
     401             : 
     402             :   // Initialize model data to zero, so corruption contains
     403             :   //  only the AMueller solution
     404             :   //  TBD: may wish to make this user togglable.
     405         136 :   vb.setModelVisCube(Complex(0.0));
     406             : 
     407         136 :   if(fitorder_p == 0){
     408             :     // Call general version:
     409         102 :     VisMueller::corrupt(vb);
     410             :   }
     411             :   else{
     412             :     // Ensure weight calibration off internally for corrupt
     413             :     //   (corruption doesn't re-scale the data!)
     414          34 :     Bool userCalWt=calWt();
     415          34 :     calWt()=false;
     416             : 
     417             :     // Bring calibration up-to-date with the vb, with inversion turned OFF
     418          34 :     syncCal(vb,false);
     419             : 
     420          34 :     Int cspw = currSpw();
     421          34 :     VBContinuumSubtractor vbcs;
     422          68 :     if (lofreq_p[cspw]>0.0 &&
     423          34 :         hifreq_p[cspw]>0.0)
     424          34 :       vbcs.init(currCPar().shape(), nAnt() - 1, totnumchan_p[cspw],
     425          34 :                 lofreq_p[cspw], hifreq_p[cspw]);
     426             :     else
     427           0 :       throw(AipsError("AMueller::applyCal: Bad freq domain info."));
     428             :     
     429          34 :     MS::PredefinedColumns whichcol = MS::MODEL_DATA;
     430          34 :     if(!vbcs.apply(vb, whichcol, currCPar(), currParOK(), false,
     431          34 :                    !spwApplied_p[cspw]))
     432             :       throw(AipsError("Could not place the continuum estimate in "
     433           0 :                       + MS::columnName(whichcol)));
     434          34 :     spwApplied_p[cspw] = true;
     435             :     // Restore user's calWt()
     436          34 :     calWt()=userCalWt; 
     437             :   }
     438         136 : }
     439             : 
     440          16 : void ANoise::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim)
     441             : {
     442          16 :   if (prtlev()>2) cout << " AN::createCorruptor()" << endl;
     443          16 :   AlwaysAssert((isSimulated()),AipsError);
     444             : 
     445          16 :   acorruptor_p = new ANoiseCorruptor();
     446          16 :   corruptor_p = acorruptor_p;
     447             : 
     448             :   // call generic parent to set corr,spw,etc info
     449          16 :   SolvableVisCal::createCorruptor(vi,simpar,nSim);
     450             : 
     451          16 :   Int Seed(1234);
     452          16 :   if (simpar.isDefined("seed")) {    
     453          16 :     Seed=simpar.asInt("seed");
     454             :   }
     455             : 
     456          16 :   Float Amp(1.0);
     457          16 :   if (simpar.isDefined("amplitude")) {    
     458          16 :     Amp=simpar.asFloat("amplitude");
     459             :   }
     460             : 
     461          16 :   acorruptor_p->initialize(Seed,Amp);
     462             : 
     463          32 :   String Mode("calc"); // calc means multiply by 1/sqrt(dnu dt)
     464          16 :   if (simpar.isDefined("mode")) {    
     465          16 :     Mode=simpar.asString("mode");
     466             :   }
     467             : 
     468          16 :   acorruptor_p->mode()=Mode;
     469             : 
     470          16 :   if (prtlev()>2) cout << " ~AN::createCorruptor()" << endl;
     471             : 
     472          16 : }
     473             : 
     474             : 
     475          16 : ANoise::ANoise(VisSet& vs) :
     476             :   VisCal(vs),             // virtual base
     477             :   VisMueller(vs),         // virtual base
     478          16 :   SolvableVisMueller(vs)  // immediate parent
     479             : {
     480          16 :   if (prtlev()>2) cout << "ANoise::ANoise(vs)" << endl;
     481          16 : }
     482             : 
     483           0 : ANoise::ANoise(String msname,Int MSnAnt,Int MSnSpw) :
     484             :   VisCal(msname,MSnAnt,MSnSpw),             // virtual base
     485             :   VisMueller(msname,MSnAnt,MSnSpw),         // virtual base
     486           0 :   SolvableVisMueller(msname,MSnAnt,MSnSpw)  // immediate parent
     487             : {
     488           0 :   if (prtlev()>2) cout << "ANoise::ANoise(msname,MSnAnt,MSnSpw)" << endl;
     489           0 : }
     490             : 
     491           0 : ANoise::ANoise(const MSMetaInfoForCal& msmc) :
     492             :   VisCal(msmc),             // virtual base
     493             :   VisMueller(msmc),         // virtual base
     494           0 :   SolvableVisMueller(msmc)  // immediate parent
     495             : {
     496           0 :   if (prtlev()>2) cout << "ANoise::ANoise(msmc)" << endl;
     497           0 : }
     498             : 
     499           0 : ANoise::ANoise(const Int& nAnt) :
     500             :   VisCal(nAnt),
     501             :   VisMueller(nAnt),
     502           0 :   SolvableVisMueller(nAnt)
     503             : {
     504           0 :   if (prtlev()>2) cout << "ANoise::ANoise(nAnt)" << endl;
     505           0 : }
     506             : 
     507          32 : ANoise::~ANoise() {
     508          16 :   if (prtlev()>2) cout << "ANoise::~ANoise()" << endl;
     509          32 : }
     510             : 
     511             : 
     512             : // Calculate Mueller matrices for all baselines
     513        4532 : void ANoise::calcAllMueller() {
     514             : 
     515        4532 :   if (prtlev()>6) cout << "        AN::calcAllMueller" << endl;
     516             : 
     517             :   // Should handle OK flags in this method, and only
     518             :   //  do Mueller calc if OK
     519             : 
     520        9064 :   Vector<Complex> oneMueller;
     521        9064 :   Vector<Bool> oneMOK;
     522        9064 :   Vector<Complex> onePar;
     523        9064 :   Vector<Bool> onePOK;
     524             : 
     525        9064 :   ArrayIterator<Complex> Miter(currMElem(),1);
     526        9064 :   ArrayIterator<Bool>    MOKiter(currMElemOK(),1);
     527        9064 :   ArrayIterator<Complex> Piter(currCPar(),1);
     528        9064 :   ArrayIterator<Bool>    POKiter(currParOK(),1);
     529             :   
     530             :   // All required baselines
     531        4532 :   Int iant1(0), iant2(-1);
     532       19082 :   for (Int ibln=0; ibln<nCalMat(); ibln++) {
     533             :     // update antenna numbers of the current baseline (assumes iant1 <= iant2)
     534       14550 :     iant2++;
     535       14550 :     if (iant2 == nAnt()) {
     536         543 :       iant1++;
     537         543 :       iant2 = iant1;
     538             :     }
     539       14550 :     corruptor_p->currAnt()=iant1;
     540       14550 :     corruptor_p->currAnt2()=iant2;
     541             : 
     542             :     // The following assumes that nChanPar()=1 or nChanMat()
     543             : 
     544       29100 :     for (Int ich=0; ich<nChanMat(); ich++) {
     545             : 
     546       14550 :       oneMueller.reference(Miter.array());
     547       14550 :       oneMOK.reference(MOKiter.array());
     548       14550 :       onePar.reference(Piter.array());
     549       14550 :       onePOK.reference(POKiter.array());
     550             :       
     551             :       // TBD  What if calcOneMueller needs freq value info?
     552             :       
     553       14550 :       calcOneMueller(oneMueller,oneMOK,onePar,onePOK);
     554             : 
     555             :       // Advance iterators, as required
     556       14550 :       Miter.next();
     557       14550 :       MOKiter.next();
     558       14550 :       if (freqDepPar()) {
     559       14550 :         Piter.next();
     560       14550 :         POKiter.next();
     561             :       }
     562             : 
     563             :     }
     564             : 
     565             :     // Step to next baseline's pars if we didn't in channel loop
     566       14550 :     if (!freqDepPar()) {
     567           0 :       Piter.next();
     568           0 :       POKiter.next();
     569             :     }
     570             :   }
     571             : 
     572        4532 : }
     573             : 
     574             : 
     575             : 
     576             : 
     577       14550 : void ANoise::calcOneMueller(Vector<Complex>& mat, Vector<Bool>& mOk,
     578             :                             const Vector<Complex>& par, const Vector<Bool>& pOk) {
     579             :   
     580       14550 :   if(prtlev() > 10)
     581             :     cout << "        AN::calcOneMueller()\n"
     582           0 :          << "par: " << par                     // These are more to remove a compiler
     583           0 :          << "\npOk: " << pOk                   // warning (par & pOK unused)
     584           0 :          << endl;                              // than anything.
     585             : 
     586             :   // If Mueller matrix is trivial, shouldn't get here
     587       14550 :   if (trivialMuellerElem()) 
     588           0 :     throw(AipsError("Trivial Mueller Matrix logic error."));
     589             :   else {
     590       14550 :     Int len=0;
     591       14550 :     mat.shape(len);
     592       43650 :     for (Int i=0; i<len; i++) {
     593       29100 :       mat[i]=acorruptor_p->simPar(); // single complex #
     594       29100 :       mOk[i]=true;
     595             :     }    
     596             :   }
     597       14550 : }
     598             : 
     599             : 
     600             : 
     601             : } //# NAMESPACE CASA - END

Generated by: LCOV version 1.16