Line data    Source code

       1             : // -*- C++ -*-
       2             : //# EVLAAperture.cc: Implementation of the EVLAAperture 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: aips2-request@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/Utils.h>
      30             : #include <synthesis/TransformMachines/EVLAAperture.h>
      31             : #include <synthesis/TransformMachines/SynthesisError.h>
      32             : #include <synthesis/TransformMachines/BeamCalc.h>
      33             : #include <synthesis/TransformMachines/WTerm.h>
      34             : #include <synthesis/TransformMachines/VLACalcIlluminationConvFunc.h>
      35             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      36             : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
      37             : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
      38             : //
      39             : //---------------------------------------------------------------------
      40             : //---------------------------------------------------------------------
      41             : // TEMPS
      42             : //---------------------------------------------------------------------
      43             : //---------------------------------------------------------------------
      44             : //
      45             : // #define CONVSIZE (1024*2)
      46             : // #define CONVWTSIZEFACTOR 1.0
      47             : // #define OVERSAMPLING 20
      48             : // #define THRESHOLD 1E-4
      49             : 
      50             : using namespace casacore;
      51             : namespace casa{
      52             :   
      53           0 :   EVLAAperture& EVLAAperture::operator=(const EVLAAperture& other)
      54             :   {
      55           0 :     if(this!=&other) 
      56             :       {
      57             :         //      ConvolutionFunction::operator=(other);
      58           0 :         logIO_p = other.logIO_p;
      59             :         //      setParams(other.polMap_p_base, other.feedStokes_p);
      60           0 :         setPolMap(other.polMap_p_base);
      61           0 :         Diameter_p=other.Diameter_p;
      62           0 :         Nant_p=other.Nant_p;
      63           0 :         HPBW=other.HPBW;
      64           0 :         sigma=other.sigma;
      65             :       }
      66           0 :     return *this;
      67             :   }
      68           0 :   Int EVLAAperture::getVLABandID(Double& vbRefFreq,String&telescopeName, const CoordinateSystem& skyCoord)
      69             :   {
      70           0 :     LogIO log_l(LogOrigin("EVLAAperture", "getVLABandID[R&D]"));
      71             : 
      72           0 :     Double refFreq = skyCoord.spectralCoordinate(skyCoord.findCoordinate(Coordinate::SPECTRAL)).referenceValue()(0);
      73           0 :     cerr << "getVLABand (Global VB Ref. min, CF Ref.): " << vbRefFreq << " , " << refFreq << endl;
      74           0 :     if (telescopeName=="VLA")
      75             :       {
      76           0 :         if      ((refFreq >=1.34E9) && (refFreq <=1.73E9)) return BeamCalc_VLA_L;
      77           0 :         else if ((refFreq >=4.5E9)  && (refFreq <=5.0E9))  return BeamCalc_VLA_C;
      78           0 :         else if ((refFreq >=8.0E9)  && (refFreq <=8.8E9))  return BeamCalc_VLA_X;
      79           0 :         else if ((refFreq >=14.4E9) && (refFreq <=15.4E9)) return BeamCalc_VLA_U;
      80           0 :         else if ((refFreq >=22.0E9) && (refFreq <=24.0E9)) return BeamCalc_VLA_K;
      81           0 :         else if ((refFreq >=40.0E9) && (refFreq <=50.0E9)) return BeamCalc_VLA_Q;
      82           0 :         else if ((refFreq >=100E6)  && (refFreq <=300E6))  return BeamCalc_VLA_4;
      83             :       }
      84             :     else 
      85           0 :       if (telescopeName=="EVLA")
      86             :         {
      87           0 :           if      ((refFreq >= 0.9E9) && (refFreq <= 2.1E9))  return BeamCalc_EVLA_L;
      88           0 :           else if ((refFreq >=2.0E9)  && (refFreq <=4.0E9))   return BeamCalc_EVLA_S;
      89           0 :           else if ((refFreq >=4.0E9)  && (refFreq <=8.0E9))   return BeamCalc_EVLA_C;
      90           0 :           else if ((refFreq >=8.0E9)  && (refFreq <=12.0E9))  return BeamCalc_EVLA_X;
      91           0 :           else if ((refFreq >=12.0E9) && (refFreq <=18.0E9))  return BeamCalc_EVLA_U;
      92           0 :           else if ((refFreq >=18.0E9) && (refFreq <=26.5E9))  return BeamCalc_EVLA_K;
      93           0 :           else if ((refFreq >=26.5E9) && (refFreq <=40.8E9))  return BeamCalc_EVLA_A;
      94           0 :           else if ((refFreq >=40.0E9) && (refFreq <=50.0E9))  return BeamCalc_EVLA_Q;
      95             :         }
      96           0 :     ostringstream mesg;
      97           0 :     log_l << telescopeName << "/" << refFreq << "(Hz) combination not recognized." << LogIO::EXCEPTION;
      98           0 :     return -1;
      99             :   }
     100             :   
     101           0 :   void EVLAAperture::cacheVBInfo(const String& telescopeName, 
     102             :                                  const Float& diameter)
     103             :   {
     104           0 :     telescopeName_p=telescopeName;
     105           0 :     Diameter_p=diameter;
     106           0 :   }
     107             : 
     108           0 :   void EVLAAperture::cacheVBInfo(const VisBuffer& vb)
     109             :   {
     110           0 :     Vector<String> telescopeNames=vb.msColumns().observation().telescopeName().getColumn();
     111           0 :     for(uInt nt=0;nt<telescopeNames.nelements();nt++)
     112             :       {
     113           0 :         if ((telescopeNames(nt) != "VLA") && (telescopeNames(nt) != "EVLA"))
     114             :           {
     115           0 :             String mesg="We can handle only (E)VLA antennas for now.\n";
     116           0 :             mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
     117           0 :             SynthesisError err(mesg);
     118           0 :             throw(err);
     119             :           }
     120           0 :         if (telescopeNames(nt) != telescopeNames(0))
     121             :           {
     122           0 :             String mesg="We do not (yet) handle multiple telescopes for A-Projection!\n";
     123           0 :             mesg += "Not yet a \"priority\"!!";
     124           0 :             SynthesisError err(mesg);
     125           0 :             throw(err);
     126             :           }
     127             :       }
     128           0 :     telescopeName_p=telescopeNames[0];
     129             : 
     130             :     //    MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
     131             :     //Freq = vb.msColumns().spectralWindow().refFrequency()(0);
     132           0 :     Diameter_p=0;
     133           0 :     Nant_p     = vb.msColumns().antenna().nrow();
     134           0 :     for (Int i=0; i < Nant_p; i++)
     135           0 :       if (!vb.msColumns().antenna().flagRow()(i))
     136             :         {
     137           0 :           Diameter_p = vb.msColumns().antenna().dishDiameter()(i);
     138           0 :           break;
     139             :         }
     140           0 :     if (Diameter_p == 0)
     141             :       {
     142           0 :         logIO() << LogOrigin("EVLAAperture", "cacheVBInfo")
     143             :                 << "No valid or finite sized antenna found in the antenna table. "
     144             :                 << "Assuming diameter = 25m."
     145             :                 << LogIO::WARN
     146           0 :                 << LogIO::POST;
     147           0 :         Diameter_p=25.0;
     148             :       }
     149           0 :     cacheVBInfo(telescopeNames[0], Diameter_p);
     150           0 :   }
     151             : 
     152           0 :   Int EVLAAperture::getBandID(const Double& freq, const String& telescopeName,const String& bandName)
     153             :   {
     154           0 :     Int bandID=0;
     155           0 :     if (!isNoOp())
     156           0 :       bandID = BeamCalc::Instance()->getBandID(freq,telescopeName,bandName);
     157             :     
     158           0 :     return bandID;
     159             :   };
     160             : 
     161           0 :   int EVLAAperture::getVisParams(const VisBuffer& vb,const CoordinateSystem& /*im*/)
     162             :   {
     163             :     Double Freq;
     164           0 :     cacheVBInfo(vb);
     165             :     
     166           0 :     Freq = vb.msColumns().spectralWindow().refFrequency()(0);
     167           0 :     Double Lambda=C::c/Freq;
     168           0 :     HPBW = Lambda/(Diameter_p*sqrt(log(2.0)));
     169           0 :     sigma = 1.0/(HPBW*HPBW);
     170             :     //    awEij.setSigma(sigma);
     171             :     //    Int bandID = getVLABandID(Freq,telescopeNames(0),im);
     172           0 :     return getBandID(Freq, telescopeName_p,"");
     173             :     // Int bandID=0;
     174             :     // if (!isNoOp())
     175             :     //   bandID = BeamCalc::Instance()->getBandID(Freq,telescopeName_p);
     176             :     
     177             :     // return bandID;
     178             :   }
     179             :   
     180           0 :   Int EVLAAperture::makePBPolnCoords(const VisBuffer&vb,
     181             :                                      const Int& convSize,
     182             :                                      const Int& convSampling,
     183             :                                      const CoordinateSystem& skyCoord,
     184             :                                      const Int& skyNx, const Int& /*skyNy*/,
     185             :                                      CoordinateSystem& feedCoord)
     186             :   //                                 Vector<Int>& cfStokes)
     187             :   {
     188           0 :     feedCoord = skyCoord;
     189             :     //
     190             :     // Make a two dimensional image to calculate auto-correlation of
     191             :     // the ideal illumination pattern. We want this on a fine grid in
     192             :     // the UV plane
     193             :     //
     194           0 :     Int directionIndex=skyCoord.findCoordinate(Coordinate::DIRECTION);
     195             : //    cout<<"Direction index is "<< directionIndex<<endl;
     196           0 :     AlwaysAssert(directionIndex>=0, AipsError);
     197           0 :     DirectionCoordinate dc=skyCoord.directionCoordinate(directionIndex);
     198           0 :     Vector<Double> sampling;
     199           0 :     sampling = dc.increment();
     200             : //    cout<<"Image sampling set to : "<<sampling<<endl;
     201           0 :     sampling*=Double(convSampling);
     202           0 :     sampling*=Double(skyNx)/Double(convSize);
     203           0 :     dc.setIncrement(sampling);
     204             : //    cout<<"Resized sampling is : "<<sampling<<endl;
     205             :     
     206           0 :     Vector<Double> unitVec(2);
     207           0 :     unitVec=convSize/2;
     208           0 :     dc.setReferencePixel(unitVec);
     209             :     
     210             :     // Set the reference value to that of the image
     211           0 :     feedCoord.replaceCoordinate(dc, directionIndex);
     212             : 
     213             :     //
     214             :     // Make an image with circular polarization axis.
     215             :     //
     216           0 :     Int NPol=0,M,N=0;
     217           0 :     M=polMap_p_base.nelements();
     218           0 :     for(Int i=0;i<M;i++) if (polMap_p_base(i) > -1) NPol++;
     219           0 :     Vector<Int> poln(NPol);
     220             :     
     221             :     Int index;
     222           0 :     Vector<Int> inStokes;
     223           0 :     index = feedCoord.findCoordinate(Coordinate::STOKES);
     224           0 :     inStokes = feedCoord.stokesCoordinate(index).stokes();
     225           0 :     N = 0;
     226             :     try
     227             :       {
     228             :         //      cerr << "### " << polMap_p_base << " " << vb.corrType() << endl;
     229           0 :         for(Int i=0;i<M;i++) if (polMap_p_base(i) > -1) {poln(N) = vb.corrType()(i);N++;}
     230           0 :         StokesCoordinate polnCoord(poln);
     231           0 :         Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
     232           0 :         feedCoord.replaceCoordinate(polnCoord,StokesIndex);
     233             :         //      cfStokes = poln;
     234             :       }
     235           0 :     catch(AipsError& x)
     236             :       {
     237             :         throw(SynthesisFTMachineError("Likely cause: Discrepancy between the poln. "
     238           0 :                                       "axis of the data and the image specifications."));
     239             :       }
     240             :     
     241           0 :     return NPol;
     242             :   }
     243             :   
     244           0 :   Bool EVLAAperture::findSupport(Array<Complex>& func, Float& threshold,Int& origin, Int& R)
     245             :   {
     246             :     Double NSteps;
     247           0 :     Int PixInc=1;
     248           0 :     Vector<Complex> vals;
     249           0 :     IPosition ndx(4,origin,0,0,0);
     250           0 :     Bool found=false;
     251           0 :     IPosition cfShape=func.shape();
     252           0 :     Int convSize = cfShape(0);
     253           0 :     for(R=convSize/4;R>1;R--)
     254             :       {
     255           0 :         NSteps = 90*R/PixInc; //Check every PixInc pixel along a
     256             :         //circle of radious R
     257           0 :         vals.resize((Int)(NSteps+0.5));
     258           0 :         vals=0;
     259           0 :         for(Int th=0;th<NSteps;th++)
     260             :           {
     261           0 :             ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
     262           0 :             ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
     263             :             
     264           0 :             if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
     265           0 :               vals(th)=func(ndx);
     266             :           }
     267           0 :         if (max(abs(vals)) > threshold)
     268           0 :           {found=true;break;}
     269             :       }
     270           0 :     return found;
     271             :   }
     272             :   
     273           0 :   void EVLAAperture::makeFullJones(ImageInterface<Complex>& pbImage,
     274             :                                    const VisBuffer& vb,
     275             :                                    Bool doSquint, Int& bandID, Double freqVal)
     276             :   {
     277             : 
     278           0 :     if (!isNoOp())
     279             :       {
     280           0 :         VLACalcIlluminationConvFunc vlaPB;
     281           0 :         Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
     282           0 :         vlaPB.setMaximumCacheSize(cachesize);
     283           0 :         bandID = getBandID(freqVal,telescopeName_p,"");
     284             :         //bandID=getVisParams(vb,pbImage.coordinates());
     285           0 :         vlaPB.makeFullJones(pbImage,vb, doSquint, bandID, freqVal);
     286             :       }
     287           0 :   }
     288             : 
     289           0 :   void EVLAAperture::applySky(ImageInterface<Complex>& outImages,
     290             :                               const VisBuffer& vb, 
     291             :                               const Bool doSquint,
     292             :                               const Int& cfKey,
     293             :                               const Int& muellerTerm,
     294             :                               const Double freqVal)
     295             :   {
     296             : //     (void)cfKey;
     297             : //     if (!isNoOp())
     298             : //       {
     299             : //      VLACalcIlluminationConvFunc vlaPB;
     300             : //      Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
     301             : //      vlaPB.setMaximumCacheSize(cachesize);
     302             : //      Int bandID;//=getVisParams(vb,outImages.coordinates());
     303             : //      bandID = getBandID(freqVal,telescopeName_p);
     304             : // //   cout<<"EVLAAperture : muellerTerm"<<muellerTerm <<"\n";
     305             : //      //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
     306             : //      Double pa=getPA(vb);
     307             : //      vlaPB.applyPB(outImages, pa, doSquint,bandID,muellerTerm,freqVal);
     308             :         
     309             : //       }
     310           0 :     Double pa=getPA(vb);
     311           0 :     applySky(outImages,pa,doSquint,cfKey,muellerTerm, freqVal);
     312           0 :   }
     313             :   //
     314             :   // This one without VB.  Should become the default (and the only one!).
     315             :   //
     316           0 :   void EVLAAperture::applySky(ImageInterface<Complex>& outImages,
     317             :                               const Double& pa,
     318             :                               const Bool doSquint,
     319             :                               const Int& cfKey,
     320             :                               const Int& muellerTerm,
     321             :                               const Double freqVal)
     322             :   {
     323             :     (void)cfKey;
     324           0 :     if (!isNoOp())
     325             :       {
     326           0 :         VLACalcIlluminationConvFunc vlaPB;
     327           0 :         Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
     328           0 :         vlaPB.setMaximumCacheSize(cachesize);
     329             :         Int bandID;//=getVisParams(vb,outImages.coordinates());
     330             :         //cout<<"EVLAAperture : muellerTerm"<<muellerTerm <<" " << telescopeName_p << endl;
     331           0 :         bandID = getBandID(freqVal,telescopeName_p,"");
     332             :         //vlaPB.applyPB(outImages, doSquint,bandID,muellerTerm,freqVal);
     333           0 :         Double pa_l=pa;  // Due to goofup in making sure complier type checking does not come in the way!
     334           0 :         vlaPB.applyPB(outImages, pa_l, doSquint,bandID,muellerTerm,freqVal);
     335             :         
     336             :       }
     337           0 :   }
     338             : 
     339           0 :   void EVLAAperture::applySky(ImageInterface<Float>& outImages,
     340             :                               const VisBuffer& vb, 
     341             :                               const Bool doSquint,
     342             :                               const Int& cfKey,
     343             :                               const Int& muellerTerm,
     344             :                               const Double freqVal)
     345             :   {
     346             :     (void)cfKey;
     347             :     (void)muellerTerm;
     348           0 :     if (!isNoOp())
     349             :       {
     350           0 :         VLACalcIlluminationConvFunc vlaPB;
     351           0 :         Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
     352           0 :         vlaPB.setMaximumCacheSize(cachesize);
     353             :         Int bandID;//=getVisParams(vb,outImages.coordinates());
     354           0 :         bandID = getBandID(freqVal,telescopeName_p,"");
     355           0 :         Double pa=getPA(vb);
     356           0 :         vlaPB.applyPB(outImages, pa, bandID, doSquint,freqVal);
     357             :       }
     358           0 :   }
     359             :   
     360             : };