Line data    Source code

       1             : // -*- C++ -*-
       2             : //# PolOuterProduct.cc: Implementation of the PolOuterProduct class
       3             : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
       4             : //# Associated Universities, Inc. Washington DC, USA.
       5             : //#
       6             : //# This library is free software; you can redistribute it and/or modify it
       7             : //# under the terms of the GNU Library General Public License as published by
       8             : //# the Free Software Foundation; either version 2 of the License, or (at your
       9             : //# option) any later version.
      10             : //#
      11             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      12             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      13             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      14             : //# License for more details.
      15             : //#
      16             : //# You should have received a copy of the GNU Library General Public License
      17             : //# along with this library; if not, write to the Free Software Foundation,
      18             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      19             : //#
      20             : //# Correspondence concerning AIPS++ should be addressed as follows:
      21             : //#        Internet email: casa-feedback@nrao.edu.
      22             : //#        Postal address: AIPS++ Project Office
      23             : //#                        National Radio Astronomy Observatory
      24             : //#                        520 Edgemont Road
      25             : //#                        Charlottesville, VA 22903-2475 USA
      26             : //#
      27             : //# $Id$
      28             : //
      29             : #include <synthesis/TransformMachines/PolOuterProduct.h>
      30             : 
      31             : using namespace casacore;
      32             : namespace casa{
      33             :   //---------------------------------------------------------------------
      34             :   //
      35             :   // Map indices to the vis. pols. used.  This remains fixed for a
      36             :   // given imaging run.
      37             :   //
      38             :   // cfIndices_p are filled by this method and is determined by the
      39             :   // Stokes-I setting for a given imaging run.  Since this also
      40             :   // depends on the Pol. mapping in the MS, this mapping is determined
      41             :   // once (but has to be determined programmatically for each MS and
      42             :   // can't be a hard-coded convention here).
      43             :   //
      44         324 :   void PolOuterProduct::initCFMaps(const Vector<Int>& vbPol, const Vector<Int>& visPol2ImMap)
      45             :   {
      46         324 :     Int np=0;
      47             :     //muellerType_p=HYBRID;
      48         942 :     for (uInt i=0;i<visPol2ImMap.nelements(); i++) 
      49             :       {
      50             :         Int n;
      51         618 :         if (visPol2ImMap(i) >= 0) 
      52             :           {
      53         612 :             n=translateStokesToGeneric(vbPol(i));
      54         612 :             cfIndices_p.resize(n+1,true);
      55         612 :             if (muellerType_p == DIAGONAL)        cfIndices_p(n).resize(1);
      56           0 :             else if (muellerType_p == FULL)       cfIndices_p(n).resize(4);
      57             :         
      58        1224 :             for(uInt j=0;j<cfIndices_p(n).nelements();j++) cfIndices_p(n)(j)=np++;
      59             :                 
      60             :           }
      61             :       }   
      62             :     //  if (visPol2ImMap(i) >= 0)
      63             :         {
      64         324 :                 if(muellerType_p == HYBRID)
      65             :                 {
      66             :                   //Int n;
      67             :                   //n=translateStokesToGeneric(vbPol(i));
      68             :                   //cerr<<"N = :"<<n<<endl;
      69           0 :                   cfIndices_p.resize(4);
      70           0 :                   cfIndices_p(0).resize(3);
      71           0 :                   cfIndices_p(1).resize(3);
      72           0 :                   cfIndices_p(2).resize(3);
      73           0 :                   cfIndices_p(3).resize(3);
      74             :                   //              cfIndices_p(0)(0) = 0;
      75             :                   //              cfIndices_p(0)(1) = 1;
      76             :                   //              cfIndices_p(0)(2) = 2;
      77             :                   //              cfIndices_p(1)(0) = 4;
      78             :                   //              cfIndices_p(1)(1) = 5;
      79             :                   //              cfIndices_p(1)(2) = 7;
      80             :                   //              cfIndices_p(2)(0) = 8;
      81             :                   //              cfIndices_p(2)(1) = 10;
      82             :                   //              cfIndices_p(2)(2) = 11;
      83             :                   //              cfIndices_p(3)(0) = 13;
      84             :                   //              cfIndices_p(3)(1) = 14;
      85             :                   //              cfIndices_p(3)(2) = 15;
      86             : 
      87           0 :                   cfIndices_p(0)(0) = 0;
      88           0 :                   cfIndices_p(0)(1) = 1;
      89           0 :                   cfIndices_p(0)(2) = 2;
      90           0 :                   cfIndices_p(1)(0) = 3;
      91           0 :                   cfIndices_p(1)(1) = 4;
      92           0 :                   cfIndices_p(1)(2) = 5;
      93           0 :                   cfIndices_p(2)(0) = 6;
      94           0 :                   cfIndices_p(2)(1) = 7;
      95           0 :                   cfIndices_p(2)(2) = 8;
      96           0 :                   cfIndices_p(3)(0) = 9;
      97           0 :                   cfIndices_p(3)(1) = 10;
      98           0 :                   cfIndices_p(3)(2) = 11;
      99             :                 }
     100             : 
     101             :         }
     102             :         //cout<<"cfIndices_p in PolOuterProduct::initCFMaps"<<cfIndices_p;
     103             : 
     104         324 :   }
     105             :   //---------------------------------------------------------------------
     106             :   //
     107             :   // Make a map between the vbPol indices and indices to get the CFs
     108             :   // to be used.  This can be dynamic, depending on the meaning of the
     109             :   // pol. axis of the visCube in the VB (which can change per VB).
     110             :   //
     111           0 :   PolMapType& PolOuterProduct::makePol2CFMat_p(const Vector<Int>& vbPol, 
     112             :                                                const Vector<Int>& vbPol2ImMap,
     113             :                                                PolMapType& outerProdNdx2VBPolMap)
     114             :   {
     115           0 :     Int nVBPol=vbPol.nelements();
     116           0 :     outerProdNdx2VBPolMap.resize(nVBPol);
     117           0 :     for (Int i=0;i<nVBPol; i++)
     118           0 :       if (vbPol2ImMap(i) >= 0)
     119             :         {
     120           0 :           Int n = translateStokesToGeneric(vbPol(i));
     121           0 :           outerProdNdx2VBPolMap(i).assign(cfIndices_p(n));
     122             :           //cout<<"cfIndices_p"<<cfIndices_p<<endl;
     123           0 :           if (muellerType_p == HYBRID)
     124             :             {
     125           0 :                 outerProdNdx2VBPolMap(i).resize(3);
     126           0 :                 Int n=translateStokesToGeneric(vbPol(i));
     127           0 :                 if (n==GPP)
     128             :                 {
     129           0 :                         outerProdNdx2VBPolMap(0)(0)=0; //cfIndices_p(0)(0)=0;
     130           0 :                         outerProdNdx2VBPolMap(0)(1)=1; //cfIndices_p(0)(1)=1;
     131           0 :                         outerProdNdx2VBPolMap(0)(2)=2; //cfIndices_p(0)(2)=2;
     132             :                         //cout<<"\n outerProdNdx2VBPolMap: "<< outerProdNdx2VBPolMap(0);
     133             :                         //cout<<"\n cfIndices_p: "<< cfIndices_p(0)<< "\n";
     134             :                 }
     135           0 :                 else if (n==GPQ)
     136             :                 {
     137           0 :                         outerProdNdx2VBPolMap(1)(0)=4; //cfIndices_p(1)(0)=4; 
     138           0 :                         outerProdNdx2VBPolMap(1)(1)=5; //cfIndices_p(1)(1)=5; 
     139           0 :                         outerProdNdx2VBPolMap(1)(2)=7; //cfIndices_p(1)(2)=7;
     140             : 
     141             :                         //cout<<"\n outerProdNdx2VBPolMap: "<< outerProdNdx2VBPolMap(1);
     142             :                         //cout<<"\n cfIndices_p: "<< cfIndices_p<< "\n";
     143             :                 }
     144           0 :                 else if (n==GQP)
     145             :                 {
     146           0 :                         outerProdNdx2VBPolMap(2)(0)=8; //cfIndices_p(2)(0)=8;
     147           0 :                         outerProdNdx2VBPolMap(2)(1)=10; //cfIndices_p(2)(1)=10;
     148           0 :                         outerProdNdx2VBPolMap(2)(2)=11; //cfIndices_p(2)(2)=11;
     149             :                         //cout<<"\n outerProdNdx2VBPolMap: "<< outerProdNdx2VBPolMap(2);
     150             :                         //cout<<"\n cfIndices_p: "<< cfIndices_p(2)<< "\n";
     151             :                 }
     152           0 :                 else if (n==GQQ)
     153             :                 {
     154           0 :                         outerProdNdx2VBPolMap(3)(0)=13; //cfIndices_p(3)(0)=13; 
     155           0 :                         outerProdNdx2VBPolMap(3)(1)=14; //cfIndices_p(3)(1)=14;
     156           0 :                         outerProdNdx2VBPolMap(3)(2)=15; //cfIndices_p(3)(2)=15;
     157             :                         //cout<<"\n outerProdNdx2VBPolMap: "<< outerProdNdx2VBPolMap(3);
     158             :                         //cout<<"\n cfIndices_p: "<< cfIndices_p(3)<< "\n";
     159             :                 }
     160             :             }
     161             :         }
     162           0 :     return outerProdNdx2VBPolMap;
     163             :   }
     164             : 
     165           0 :   PolMapType& PolOuterProduct::makePol2CFMat(const Vector<Int>& vbPol, 
     166             :                                              const Vector<Int>& vbPol2ImMap)
     167           0 :   {return makePol2CFMat_p(vbPol,vbPol2ImMap, outerProductIndex2VBPolMap_p);}
     168             : 
     169           0 :   PolMapType& PolOuterProduct::makeConjPol2CFMat(const Vector<Int>& vbPol, 
     170             :                                                  const Vector<Int>& /*vbPol2ImMap*/)
     171             :   {
     172           0 :     Int nVBPol=vbPol.nelements();
     173             :     // Resize the conj. index map.
     174           0 :     conjOuterProductIndex2VBPolMap_p.resize(nVBPol);
     175           0 :     for (Int i=0;i<nVBPol;i++)
     176           0 :       conjOuterProductIndex2VBPolMap_p(i).assign(outerProductIndex2VBPolMap_p(i));
     177             :     //
     178             :     // For each conjOuterProduct entry, find the equivalent entry in outerProductMap.
     179             :     // Copy of the associated index from outerProductIndex map to conjOuterProductMap
     180             :     //
     181           0 :     for (uInt i=0;i<conjOuterProduct2VBPolMap_p.nelements();i++)
     182           0 :       for (uInt j=0;j<conjOuterProduct2VBPolMap_p(i).nelements();j++)
     183             :         {
     184           0 :           Int el=conjOuterProduct2VBPolMap_p(i)(j);
     185             : 
     186           0 :           for (uInt ii=0;ii<outerProduct2VBPolMap_p.nelements();ii++)
     187           0 :             for (uInt jj=0;jj<outerProduct2VBPolMap_p(ii).nelements();jj++)
     188           0 :               if (outerProduct2VBPolMap_p(ii)(jj) == el)
     189             :                 {
     190           0 :                   conjOuterProductIndex2VBPolMap_p(i)(j) = outerProductIndex2VBPolMap_p(ii)(jj); 
     191           0 :                   break;
     192             :                 }
     193             :         }
     194           0 :     return conjOuterProductIndex2VBPolMap_p;
     195             :     //    return makePol2CFMat_p(vbPol,vbPol2ImMap, conjOuterProductIndex2VBPolMap_p);
     196             :   }
     197             : 
     198             :   //---------------------------------------------------------------------
     199             :   //
     200             :   // Map the CF indices to physical outer product enums.
     201             :   //
     202             :   // The muellerRows_p mapping is an internal convention and is hard coded in protected method makeMap().
     203             :   //
     204           0 :   PolMapType& PolOuterProduct::makePolMat_p(const Vector<Int>& vbPol, 
     205             :                                             const Vector<Int>& vbPol2ImMap,
     206             :                                             PolMapType& outerProd2VBPolMap,
     207             :                                             RigidVector<RigidVector<Int, 4>, 4>& mRows)
     208             :   {
     209           0 :     Int nVBPol=vbPol.nelements();
     210             :     
     211           0 :     outerProd2VBPolMap.resize(nVBPol);
     212           0 :     for (Int i=0; i<nVBPol; i++)
     213             :       {
     214           0 :         if (vbPol2ImMap(i) >= 0)
     215             :           {
     216           0 :             if (muellerType_p == DIAGONAL)
     217             :               {
     218           0 :                 outerProd2VBPolMap(i).resize(1);
     219           0 :                 Int n=translateStokesToGeneric(vbPol(i));
     220           0 :                 if (n==GPP) outerProd2VBPolMap(i)(0)=mRows(n)(0);
     221           0 :                 if (n==GPQ) outerProd2VBPolMap(i)(0)=mRows(n)(1);
     222           0 :                 if (n==GQP) outerProd2VBPolMap(i)(0)=mRows(n)(2);
     223           0 :                 if (n==GQQ) outerProd2VBPolMap(i)(0)=mRows(n)(3);
     224             :               }
     225           0 :             else if (muellerType_p == FULL)
     226             :               {
     227           0 :                 outerProd2VBPolMap(i).resize(4);
     228           0 :                 Int n=translateStokesToGeneric(vbPol(i));
     229           0 :                 outerProd2VBPolMap(i).assign(mRows(n).vector());
     230             :         //      cout<<"\n outerProd2VBPolMap: "<< outerProd2VBPolMap(i);
     231             :         //      cout<<"\n mRows: "<< mRows(n)<< "\n";
     232             :               }
     233             : // The Hybrid option ensures that we are omitting the anti-diagonal that is posing some normalization issue
     234             : // We do this by ensuring that the outerProd2VbPolMap is resized to 3 and set with the appropriate index.
     235           0 :             else if (muellerType_p == HYBRID)
     236             :             {
     237           0 :                 outerProd2VBPolMap(i).resize(3);
     238           0 :                 Int n=translateStokesToGeneric(vbPol(i));
     239           0 :                 if (n==GPP)
     240             :                 {
     241           0 :                         outerProd2VBPolMap(0)(0)=mRows(0)(0);
     242           0 :                         outerProd2VBPolMap(0)(1)=mRows(0)(1);
     243           0 :                         outerProd2VBPolMap(0)(2)=mRows(0)(2);
     244             :                         //cout<<"\n outerProd2VBPolMap: "<< outerProd2VBPolMap(0);
     245             :                         //cout<<"\n mRows: "<< mRows(0)<< "\n";
     246             :                 }
     247           0 :                 else if (n==GPQ)
     248             :                 {
     249           0 :                         outerProd2VBPolMap(1)(0)=mRows(1)(0);
     250           0 :                         outerProd2VBPolMap(1)(1)=mRows(1)(1);
     251           0 :                         outerProd2VBPolMap(1)(2)=mRows(1)(3);
     252             :                         //cout<<"\n outerProd2VBPolMap: "<< outerProd2VBPolMap(1);
     253             :                         //cout<<"\n mRows: "<< mRows(1)<< "\n";
     254             :                 }
     255           0 :                 else if (n==GQP)
     256             :                 {
     257           0 :                         outerProd2VBPolMap(2)(0)=mRows(2)(0);
     258           0 :                         outerProd2VBPolMap(2)(1)=mRows(2)(2);
     259           0 :                         outerProd2VBPolMap(2)(2)=mRows(2)(3);
     260             :                         //cout<<"\n outerProd2VBPolMap: "<< outerProd2VBPolMap(2);
     261             :                         //cout<<"\n mRows: "<< mRows(2)<< "\n";
     262             :                 }
     263           0 :                 else if (n==GQQ)
     264             :                 {
     265           0 :                         outerProd2VBPolMap(3)(0)=mRows(3)(1);
     266           0 :                         outerProd2VBPolMap(3)(1)=mRows(3)(2);
     267           0 :                         outerProd2VBPolMap(3)(2)=mRows(3)(3);
     268             :                         //cout<<"\n outerProd2VBPolMap: "<< outerProd2VBPolMap(3);
     269             :                         //cout<<"\n mRows: "<< mRows(3)<< "\n";
     270             :                 }
     271             :             }
     272             :         
     273             :           }
     274             :       }
     275           0 :     return outerProd2VBPolMap;
     276             :   }
     277           0 :   PolMapType& PolOuterProduct::makePolMat(const Vector<Int>& vbPol, 
     278             :                                           const Vector<Int>& vbPol2ImMap)
     279           0 :   {return makePolMat_p(vbPol,vbPol2ImMap, outerProduct2VBPolMap_p, muellerRows_p);}
     280           0 :   PolMapType& PolOuterProduct::makeConjPolMat(const Vector<Int>& vbPol, 
     281             :                                           const Vector<Int>& vbPol2ImMap)
     282           0 :   {return makePolMat_p(vbPol,vbPol2ImMap, conjOuterProduct2VBPolMap_p, conjMuellerRows_p);}
     283             :   //
     284             :   //---------------------------------------------------------------------
     285             :   //
     286           0 :   void PolOuterProduct::makePolMap(const Vector<CrossPolCircular>& pols)
     287             :   {
     288           0 :     for (uInt i=0;i<pols.nelements(); i++)
     289           0 :       setElement(pols(i),i);
     290           0 :   }
     291             :   //---------------------------------------------------------------------
     292             :   //
     293             :   // ToDo:Make hardcoded maps for this
     294         612 :   PolOuterProduct::GenericVBPol translateStokesToGeneric(const Int& vbPol)
     295             :   {
     296         612 :     if ((vbPol==Stokes::RR) || (vbPol==Stokes::XX)) return PolOuterProduct::GPP;
     297         290 :     if ((vbPol==Stokes::RL) || (vbPol==Stokes::XY)) return PolOuterProduct::GPQ;
     298         290 :     if ((vbPol==Stokes::LR) || (vbPol==Stokes::YX)) return PolOuterProduct::GQP;
     299         290 :     if ((vbPol==Stokes::LL) || (vbPol==Stokes::YY)) return PolOuterProduct::GQQ;
     300           0 :     return PolOuterProduct::JUSTWRONGVBPOL;
     301             :   }
     302             :   //
     303             :   //---------------------------------------------------------------------
     304             :   // ToDo: Make hardcoded maps for this
     305           0 :   Int translateGenericToStokes(const PolOuterProduct::GenericVBPol& gPol,
     306             :                                const MSIter::PolFrame& polFrame)
     307             :   {
     308           0 :     if (polFrame==MSIter::Circular)
     309             :       {
     310           0 :         if (gPol==PolOuterProduct::GPP) return Stokes::RR;
     311           0 :         if (gPol==PolOuterProduct::GPQ) return Stokes::RL;
     312           0 :         if (gPol==PolOuterProduct::GQP) return Stokes::LR;
     313           0 :         if (gPol==PolOuterProduct::GQQ) return Stokes::LL;
     314             :       }
     315           0 :     else if (polFrame==MSIter::Linear)
     316             :       {
     317           0 :         if (gPol==PolOuterProduct::GPP) return Stokes::XX;
     318           0 :         if (gPol==PolOuterProduct::GPQ) return Stokes::XY;
     319           0 :         if (gPol==PolOuterProduct::GQP) return Stokes::YX;
     320           0 :         if (gPol==PolOuterProduct::GQQ) return Stokes::YY;
     321             :       }
     322           0 :     return PolOuterProduct::JUSTWRONGVBPOL;
     323             :   }
     324             :   
     325             : };