//# ImagePolProxy.cc:  C++ casa namespace proxy for ImagePol tool
//# Copyright (C) 2007
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//# $Id: ImagePolProxy.cc 20615 2009-06-09 02:16:01Z Malte.Marquarding $
#include <imageanalysis/ImageAnalysis/ImagePolProxy.h>
#include <casacore/casa/aips.h>
#include <iostream>
#include <sstream>
#include <casacore/casa/IO/ArrayIO.h>
#include <casacore/casa/Arrays/ArrayMath.h>
#include <casacore/casa/Arrays/ArrayUtil.h>
#include <casacore/casa/Arrays/MaskedArray.h>
#include <casacore/casa/Arrays/MaskArrMath.h>
#include <casacore/casa/BasicMath/Random.h>
#include <casacore/casa/BasicSL/String.h>
#include <casacore/casa/Containers/Record.h>
#include <casacore/casa/Exceptions/Error.h>
#include <casacore/casa/Logging/LogFilter.h>
#include <casacore/casa/Logging/LogIO.h>
#include <casacore/casa/Logging/LogOrigin.h>
#include <casacore/coordinates/Coordinates/CoordinateUtil.h>
#include <casacore/coordinates/Coordinates/StokesCoordinate.h>
#include <casacore/coordinates/Coordinates/CoordinateSystem.h>
#include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
#include <casacore/images/Images/ImageExpr.h>
#include <casacore/images/Images/ImageInterface.h>
#include <casacore/images/Regions/ImageRegion.h>
#include <casacore/images/Images/ImageUtilities.h>
#include <casacore/images/Images/PagedImage.h>
#include <casacore/images/Regions/RegionHandler.h>
#include <casacore/images/Images/SubImage.h>
#include <casacore/images/Images/TempImage.h>
#include <casacore/lattices/Lattices/LatticeUtilities.h>
#include <casacore/measures/Measures/Stokes.h>
#include <casacore/tables/Tables/Table.h>
#include <casacore/tables/LogTables/NewFile.h>
#include <casacore/casa/namespace.h>

#include <memory>

using namespace casacore;
namespace casa { //# name space casa begins


  ImagePol::ImagePol() : itsImPol(0)
  {
    itsLog= new LogIO();
    
  }

  ImagePol::ImagePol(ImageInterface<Float>& im) :  itsLog(0),itsImPol(0)
  {
    open(im);
    
  }

  ImagePol::~ImagePol(){
    if(itsLog)
      delete itsLog;
    if(itsImPol)
      delete itsImPol;
  }

  Bool ImagePol::open(ImageInterface<Float>& im){
    if(!itsLog)
      itsLog= new LogIO();
    if(itsImPol)
      delete itsImPol;
    itsImPol= new ImagePolarimetry(im);
    
    return true;

  }

  Bool ImagePol::open(const String& infile) {
    Bool rstat(false);

    if(!itsLog) itsLog = new LogIO();
    *itsLog << LogOrigin("ImagePol", "open(const String& infile)");
    if (infile.size() == 0) {
      *itsLog << LogIO::WARN << "You must give an infile" << LogIO::POST;
      return rstat;
    }

    if(itsImPol) delete itsImPol;
    itsImPol=0;
    //
    ImageInterface<Float>* imagePointer = 0;
    ImageUtilities::openImage (imagePointer, infile);

    //
    try {
      itsImPol = new ImagePolarimetry(*imagePointer);
      rstat = true;
    } catch (AipsError x) {
      delete imagePointer;
      *itsLog << x.getMesg() << LogIO::EXCEPTION;
    }
    //
    delete imagePointer;
    imagePointer = 0;
    return rstat;
  }

  Bool
  ImagePol::imagepoltestimage(const String& outFile, const Vector<Double>& rm,
			      Bool rmDefault, Double pa0, Double sigma,
			      Int nx, Int ny, Int nf, Double f0, Double df)
  {

    Bool rstat(false);

    if (itsLog==0) itsLog= new LogIO();
    *itsLog << LogOrigin("ImagePol", "imagepoltestimage");

    if (outFile.size() == 0) {
      *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;
      return rstat;
    }

    if(itsImPol) delete itsImPol;
    itsImPol = 0;

    // If not given make RM with no ambiguity
    Vector<Float> rm2(rm.size());
    if (rmDefault) {
      Double l1 = QC::c( ).getValue(Unit("m/s")) / f0;
      Double l2 = QC::c( ).getValue(Unit("m/s")) / (f0+df);
      rm2.resize(1);
      rm2(0) = C::pi / 2 / (l1*l1 - l2*l2);
    } else {
      for (uInt i = 0; i < rm.size(); i++) rm2[i] = rm[i];
    }

    const uInt nRM = rm2.nelements();
    //
    if (nRM == 1) {
      *itsLog << LogIO::NORMAL << "Using Rotation Measure = " << rm2(0)
      	      << " radians/m/m" << endl;
    } else {
      *itsLog << LogIO::NORMAL  << "Using Rotation Measures : " << endl;
      for (uInt i=0; i<nRM; i++) {
	*itsLog << "                          " << rm2(i)
		<< " radians/m/m" << endl;
      }
    }

    *itsLog << "Using pa0              = " << pa0 << " degrees" << endl;
    *itsLog << "Using frequency        = " << f0 << " Hz" << endl;
    *itsLog << "Using bandwidth        = " << df << " Hz " << endl;
    *itsLog << "Using number channels  = " << nf << LogIO::POST;

    // Make image
    IPosition shape(4,nx,ny,4,nf);
    ImageInterface<Float>* pImOut = 0;
    makeIQUVImage(pImOut, outFile, sigma, pa0*C::pi/180.0, rm2, shape, f0, df);
    try {
      itsImPol = new ImagePolarimetry(*pImOut);
      rstat = true;
    } catch (AipsError x) {
      delete pImOut;
      pImOut = 0;
      *itsLog << x.getMesg() << LogIO::EXCEPTION;
    }
    //
    delete pImOut;
    pImOut = 0;

    return rstat;
  }


  Bool ImagePol::depolratio(ImageInterface<Float>*& returnim,
			    const String& infile, Bool debias, 
			    Double clip, Double sigma, 
			    const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog <<"No attached image, please use open "
	      << LogIO::EXCEPTION;
      return rstat;
    }
    auto im1=itsImPol->imageInterface();
    std::unique_ptr<ImageInterface<Float> > im2;
    ImageUtilities::openImage(im2, infile);
    ImageExpr<Float> tmpim=itsImPol->depolarizationRatio(*im1, *im2, 
							 debias, 
							 Float(clip),
							 Float(sigma));

    String err;
    if(!copyImage(returnim, tmpim, outfile, true)){
      *itsLog <<"Could not convert ratio image "
	      << LogIO::EXCEPTION;
    }
    rstat = true;
    return rstat;
    
  }

  
Bool ImagePol::complexlinpol(const String& outfile) {
    *itsLog << LogOrigin("ImagePolProxy", __FUNCTION__);
    if (! itsImPol) {
        *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	        << LogIO::POST;
        return false;
    }

    if (outfile.size() == 0) {
        *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;
        return false;
    }
    CoordinateSystem cSysPol;
    const auto shapePol = itsImPol->singleStokesShape(
            cSysPol, Stokes::Plinear
        );
    const auto isMasked = itsImPol->isMasked();
    auto pOutComplex = _makeImage(
        outfile, cSysPol, shapePol, isMasked, false
    );
    auto expr = itsImPol->complexLinearPolarization();
    fiddleStokesCoordinate(*pOutComplex, Stokes::Plinear);

    // Copy to output

    pOutComplex->setCoordinateInfo(expr.coordinates());
    LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
    auto p = itsImPol->imageInterface();
    copyMiscellaneous (*pOutComplex, *p);
    return true;
}
  
  // Summary
  void ImagePol::summary() const {
    *itsLog << LogOrigin("imagepol", "summary()");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return;
    }
    itsImPol->summary(*itsLog);
  }

  // sigma
  Float ImagePol::sigma(Float clip) const {
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigma(clip);
  }

  // Stokes I
  Bool ImagePol::stokesI(ImageInterface<Float>*& rtnim, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "stokesI(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }
    ImageExpr<Float> expr = itsImPol->stokesI();
    // Create output image
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }

  Float ImagePol::sigmaStokesI(Float clip) const {
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaStokesI(clip);
  }

  // Stokes Q
  Bool ImagePol::stokesQ(ImageInterface<Float>*& rtnim, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "stokesQ(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }
    ImageExpr<Float> expr = itsImPol->stokesQ();

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }

  Float ImagePol::sigmaStokesQ(Float clip) const {
    *itsLog << LogOrigin("imagepol", "sigmaStokesQ(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaStokesQ(clip);
  }

  // Stokes U
  Bool ImagePol::stokesU(ImageInterface<Float>*& rtnim, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "stokesU(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }
    ImageExpr<Float> expr = itsImPol->stokesU();

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }

  Float ImagePol::sigmaStokesU(Float clip) const {
    *itsLog << LogOrigin("imagepol", "sigmaStokesU(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaStokesU(clip);
  }

  // Stokes V
  Bool ImagePol::stokesV(ImageInterface<Float>*& rtnim, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "stokesV(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }
    ImageExpr<Float> expr = itsImPol->stokesV();

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }

  Float ImagePol::sigmaStokesV(Float clip) const {
   *itsLog << LogOrigin("imagepol", "sigmaStokesV(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaStokesV(clip);
  }

  Float ImagePol::sigmaLinPolInt(Float clip, Float sigma) const {
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaLinPolInt(clip, sigma);
  }

  Float ImagePol::sigmaTotPolInt(Float clip, Float sigma) const {
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return -1.0;
    }
    return itsImPol->sigmaTotPolInt(clip, sigma);
  }

  // Complex linear polarization
  void ImagePol::complexLinearPolarization (const String& outfile) {
    *itsLog << LogOrigin("imagepol", "complexLinearPolarization(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return;
    }

    if (outfile.size() == 0)
      *itsLog << LogIO::WARN << "You must give an outfile" << LogIO::POST;

    // Make output complex image
    ImageInterface<Complex>* pOutComplex = 0;
    CoordinateSystem cSysPol;
    IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::Plinear);
    _makeImage (pOutComplex, outfile, cSysPol, shapePol,
	       itsImPol->isMasked(), false);

    // Make Expr
    ImageExpr<Complex> expr = itsImPol->complexLinearPolarization();
    fiddleStokesCoordinate(*pOutComplex, Stokes::Plinear);

    // Copy to output
    pOutComplex->setCoordinateInfo(expr.coordinates());
    LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
    //
    auto p = itsImPol->imageInterface();
    copyMiscellaneous (*pOutComplex, *p);
    delete pOutComplex;
  }

  // Complex linear polarization
  void ImagePol::complexFractionalLinearPolarization (const String& outfile) {
    *itsLog << LogOrigin("imagepol", "complexFractionalLinearPolarization(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return;
    }

    // Make output complex image
    ImageInterface<Complex>* pOutComplex = 0;
    CoordinateSystem cSysPol;
    IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::PFlinear);
    _makeImage (pOutComplex, outfile, cSysPol, shapePol,
	       itsImPol->isMasked(), false);

    std::unique_ptr<ImageInterface<Complex> > x(pOutComplex);

    // Make Expr
    ImageExpr<Complex> expr = itsImPol->complexFractionalLinearPolarization();
    fiddleStokesCoordinate(*pOutComplex, Stokes::PFlinear);

    // Copy to output
    pOutComplex->setCoordinateInfo(expr.coordinates());
    LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, expr);
    copyMiscellaneous (*pOutComplex, *(itsImPol->imageInterface()));
  }

  Bool ImagePol::sigmaLinPolPosAng(ImageInterface<Float>*& rtnim, Float clip,
				   Float sigma, const String& outfile) {
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }

    Bool radians = false;
    ImageExpr<Float> expr = itsImPol->sigmaLinPolPosAng(radians, clip, sigma);

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  // Fractional linearly polarized intensity
  Bool ImagePol::fracLinPol(ImageInterface<Float>*& rtnim, Bool debias,
			    Float clip, Float sigma, const String& outfile) {
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "fracLinPol(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }

    ImageExpr<Float> expr = itsImPol->fracLinPol(debias, clip, sigma);

    // Create output image
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  Bool ImagePol::sigmaFracLinPol(ImageInterface<Float>*& rtnim, Float clip,
				 Float sigma, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }

    ImageExpr<Float> expr = itsImPol->sigmaFracLinPol(clip, sigma);

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  // Fractional total polarized intensity
  Bool ImagePol::fracTotPol(ImageInterface<Float>*& rtnim, Bool debias,
			    Float clip, Float sigma, const String& outfile) {
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", "fracTotPol(...)");
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }

    ImageExpr<Float> expr = itsImPol->fracTotPol(debias, clip, sigma);

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  Bool ImagePol::sigmaFracTotPol(ImageInterface<Float>*& rtnim, Float clip,
		       Float sigma, const String& outfile){
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }

    ImageExpr<Float> expr = itsImPol->sigmaFracTotPol(clip, sigma);

    // Create output image if needed
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  // Depolarization ratio
  Bool ImagePol::depolarizationRatio (ImageInterface<Float>*& rtnim, 
			    const String& infile,
			    Bool debias, Float clip,
			    Float sigma, const String& outfile) {
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << "No attached image, please use open "
	      << LogIO::EXCEPTION;
    }
    auto imagePointer1 = itsImPol->imageInterface();
    ImageInterface<Float>* imagePointer2 = 0;
    ImageUtilities::openImage (imagePointer2, infile);
    //
    ImageExpr<Float> expr =
      ImagePolarimetry::depolarizationRatio(*imagePointer1, *imagePointer2,
					    debias, clip, sigma);
    //
    if (imagePointer2) delete imagePointer2;
    imagePointer2 = 0;
    
    // Create output image
    rstat = copyImage (rtnim, expr, outfile, true);

    return rstat;
  }


  Bool ImagePol::sigmaDepolarizationRatio (ImageInterface<Float>*& rtnim,
				 const String& infile,
				 Bool debias, Float clip,
				 Float sigma, const String& outfile) {
    Bool rstat(false);
    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return rstat;
    }
    //
    auto imagePointer1 = itsImPol->imageInterface();
    //
    ImageInterface<Float>* imagePointer2 = 0;
    ImageUtilities::openImage(imagePointer2, infile);
    //
    ImageExpr<Float> expr =
      ImagePolarimetry::sigmaDepolarizationRatio(*imagePointer1,
						 *imagePointer2,
						 debias, clip, sigma);
    //
    if (imagePointer2) delete imagePointer2;
    imagePointer2 = 0;

    // Create output image
    rstat = copyImage (rtnim, expr, outfile, true);
    return rstat;
  }


  // Find Rotation Measure from Fourier method
  void ImagePol::fourierRotationMeasure(const String& outfile,
			      const String& outfileAmp,
			      const String& outfilePA,
			      const String& outfileReal,
			      const String& outfileImag,
			      Bool zeroZeroLag) {

    *itsLog << LogOrigin("imagepol", __FUNCTION__, WHERE);
    if(itsImPol==0){
      *itsLog <<"No attached image, please use open "
	      << LogIO::EXCEPTION;
      return;
    }
    if (itsImPol->imageInterface()->imageInfo().hasMultipleBeams()) {
    	*itsLog << "Cannot run " << __FUNCTION__
    		<< " on an image with multiple beams"
    		<< LogIO::EXCEPTION;
    }


    // Make output complex image
    ImageInterface<Complex>* pOutComplex = 0;
    CoordinateSystem cSysPol;
    IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::Plinear);
    _makeImage (pOutComplex, outfile, cSysPol, shapePol,
	       itsImPol->isMasked(), true);

    // Make output amplitude and position angle images
    ImageInterface<Float>* pOutAmp = 0;
    ImageInterface<Float>* pOutPA = 0;
    makeImage (pOutAmp, outfileAmp, cSysPol, shapePol,
	       itsImPol->isMasked(), false);
    makeImage (pOutPA, outfilePA, cSysPol, shapePol,
	       itsImPol->isMasked(), false);

    // Make output real and imaginary images
    ImageInterface<Float>* pOutReal = 0;
    ImageInterface<Float>* pOutImag = 0;
    makeImage (pOutReal, outfileReal, cSysPol, shapePol,
	       itsImPol->isMasked(), false);
    makeImage (pOutImag, outfileImag, cSysPol, shapePol,
	       itsImPol->isMasked(), false);

    // The output complex image will have correct Coordinates, mask, and
    // miscellaneous things copied to it
    itsImPol->fourierRotationMeasure(*pOutComplex, zeroZeroLag);

    // Copy to output
    auto p = itsImPol->imageInterface();
    if (pOutAmp!=0) {
      LatticeExprNode node(abs(*pOutComplex));
      LatticeExpr<Float> le(node);
      LatticeUtilities::copyDataAndMask(*itsLog, *pOutAmp, le);
      //
      pOutAmp->setCoordinateInfo(pOutComplex->coordinates());
      copyMiscellaneous (*pOutAmp, *p);
      pOutAmp->setUnits(p->units());
      fiddleStokesCoordinate(*pOutAmp, Stokes::Plinear);
      delete pOutAmp;
    }
    if (pOutPA!=0) {
      LatticeExprNode node(pa(imag(*pOutComplex),real(*pOutComplex)));  // degrees
      LatticeExpr<Float> le(node);
      LatticeUtilities::copyDataAndMask(*itsLog, *pOutPA, le);
      //
      pOutPA->setCoordinateInfo(pOutComplex->coordinates());
      copyMiscellaneous (*pOutPA, *p);
      pOutPA->setUnits("deg");
      fiddleStokesCoordinate(*pOutPA, Stokes::Pangle);
      delete pOutPA;
    }
    if (pOutReal!=0) {
      LatticeExprNode node(real(*pOutComplex));
      LatticeExpr<Float> le(node);
      LatticeUtilities::copyDataAndMask(*itsLog, *pOutReal, le);
      pOutReal->setCoordinateInfo(pOutComplex->coordinates());
      copyMiscellaneous (*pOutReal, *p);
      pOutReal->setUnits(p->units());
      fiddleStokesCoordinate(*pOutReal, Stokes::Plinear);  // Not strictly correct
      delete pOutReal;
    }
    if (pOutImag!=0) {
      LatticeExprNode node(imag(*pOutComplex));
      LatticeExpr<Float> le(node);
      LatticeUtilities::copyDataAndMask(*itsLog, *pOutImag, le);
      pOutImag->setCoordinateInfo(pOutComplex->coordinates());
      copyMiscellaneous (*pOutImag, *p);
      pOutImag->setUnits(p->units());
      fiddleStokesCoordinate(*pOutImag, Stokes::Plinear);  // Not strictly correct
      delete pOutImag;
    }
  }

void ImagePol::rotationMeasure(
    const String& outRM, const String& outRMErr,
    const String& outPA0, const String& outPA0Err,
    const String& outNTurns, const String& outChiSq,
    Int axis2, Float sigmaQU, Float rmFg,
    Float rmMax, Float maxPaErr
) {
    // Find Rotation Measure from traditional method
    *itsLog << LogOrigin("ImagePol", __func__);
    if(itsImPol==0) {
        *itsLog << LogIO::SEVERE <<"No attached image, please use open "
            << LogIO::POST;
        return;
    }
    // Make output images.  Give them all a mask as we don't know if output
    // will be masked or not.
    CoordinateSystem cSysRM;
    Int fAxis, sAxis;
    Int axis = axis2;
    IPosition shapeRM = itsImPol->rotationMeasureShape(
        cSysRM, fAxis, sAxis, *itsLog, axis
    );
    ImageInterface<Float>* pRMOut = 0;
    ImageInterface<Float>* pRMOutErr = 0;
    makeImage (pRMOut, outRM, cSysRM, shapeRM, true, false);
    // manage naked pointers so exception throwing doesn't leave open images
    std::unique_ptr<ImageInterface<Float> > managed5(pRMOut);
    makeImage (pRMOutErr, outRMErr, cSysRM, shapeRM, true, false);
    std::unique_ptr<ImageInterface<Float> > managed6(pRMOutErr);
    CoordinateSystem cSysPA;
    IPosition shapePA = itsImPol->positionAngleShape(
        cSysPA, fAxis, sAxis, *itsLog, axis
    );
    ImageInterface<Float>* pPA0Out = 0;
    ImageInterface<Float>* pPA0OutErr = 0;
    makeImage (pPA0Out, outPA0, cSysPA, shapePA, true, false);
    std::unique_ptr<ImageInterface<Float> > managed1(pPA0Out);
    makeImage (pPA0OutErr, outPA0Err, cSysPA, shapePA, true, false);
    std::unique_ptr<ImageInterface<Float> > managed2(pPA0OutErr);
    ImageInterface<Float>* pNTurnsOut = 0;
    makeImage (pNTurnsOut, outNTurns, cSysRM, shapeRM, true, false);
    std::unique_ptr<ImageInterface<Float> > managed3(pNTurnsOut);
    ImageInterface<Float>* pChiSqOut = 0;
    makeImage (pChiSqOut, outChiSq, cSysRM, shapeRM, true, false);
    std::unique_ptr<ImageInterface<Float> > managed4(pChiSqOut);
    itsImPol->rotationMeasure(
    	pRMOut, pRMOutErr, pPA0Out, pPA0OutErr,
    	pNTurnsOut, pChiSqOut, 
    	axis, rmMax, maxPaErr,
    	sigmaQU, rmFg, true
    );
    auto p = itsImPol->imageInterface();
    if (pRMOut) {
        copyMiscellaneous (*pRMOut, *p);
    }
    if (pRMOutErr) {
        copyMiscellaneous (*pRMOutErr, *p);
    }
    if (pPA0Out) {
        copyMiscellaneous (*pPA0Out, *p);
    }
    if (pPA0OutErr) {
      copyMiscellaneous (*pPA0OutErr, *p);
    }
    if (pNTurnsOut) {
        copyMiscellaneous (*pNTurnsOut, *p);
    }
    if (pChiSqOut) {
        copyMiscellaneous (*pChiSqOut, *p);
    }
}


  // Make a complex image
  void ImagePol::makeComplex (const String& outfile, const String& real,
		    const String& imag, const String& amp,
		    const String& phase) {

    *itsLog << LogOrigin("imagepol", __FUNCTION__);
    if(itsImPol==0){
      *itsLog << LogIO::SEVERE <<"No attached image, please use open " 
	      << LogIO::POST;
      return;
    }

    // Checks
    if (outfile.empty()) {
      *itsLog << "You must give the output complex image file name"
	      << LogIO::EXCEPTION;
    }

    Bool doRI = !real.empty() && !imag.empty();
    Bool doAP = !amp.empty() && !phase.empty();
    if (doRI && doAP) {
      *itsLog << "You must give either real and imaginary, or amplitude and phase; Not all four of them"
	      << LogIO::EXCEPTION;
    }

    // Make output complex image
    ImageInterface<Complex>* pOutComplex = 0;
    CoordinateSystem cSysPol;
    IPosition shapePol = itsImPol->singleStokesShape(cSysPol, Stokes::I);
    _makeImage (pOutComplex, outfile, cSysPol, shapePol,
	       itsImPol->isMasked(), false);

    // Make Expression. Only limited Stokes types that make sense are allowed.
    ImageExpr<Complex>* pExpr = 0;
    if (doRI) {
      PagedImage<Float> rr(real);
      Stokes::StokesTypes tR = stokesType(rr.coordinates());
      //
      PagedImage<Float> ii(imag);
      Stokes::StokesTypes tI = stokesType(ii.coordinates());
      //
      if (tR!=Stokes::Q || tI!=Stokes::U) {
	*itsLog << "The real and imaginary components must be Q and U, respectively"
		<< LogIO::EXCEPTION;
      }
      Stokes::StokesTypes typeOut = Stokes::Plinear;
      //
      LatticeExprNode node(formComplex(rr,ii));
      LatticeExpr<Complex> le(node);
      pExpr = new ImageExpr<Complex>(le, String("ComplexLinearPolarization"));
      fiddleStokesCoordinate(*pExpr, typeOut);
    } else {
      PagedImage<Float> aa(amp);
      Stokes::StokesTypes tA = stokesType(aa.coordinates());
      //
      PagedImage<Float> pp(phase);
      Stokes::StokesTypes tP = stokesType(pp.coordinates());
      //
      if (tP!=Stokes::Pangle) {
	*itsLog << "The phase must be of Stokes type position angle (Pangle)"
		<< LogIO::EXCEPTION;
      }
      Float fac = 1.0;
      String units = pp.units().getName();
      if (units.contains(String("deg"))) {
	fac = C::pi / 180.0;
      } else if (units.contains(String("rad"))) {
      } else {
	*itsLog << LogIO::WARN
		<< "Units for phase are neither radians nor degrees. radians assumed"
		<< LogIO::POST;
      }
      //
      Stokes::StokesTypes typeOut = Stokes::Undefined;
      String exprName("");
      if (tA==Stokes::Ptotal) {
	typeOut = Stokes::Ptotal;
	exprName = String("ComplexTotalPolarization");
      } else if (tA==Stokes::Plinear) {
	typeOut = Stokes::Plinear;
	exprName = String("ComplexLinearPolarization");
      } else if (tA==Stokes::PFtotal) {
	typeOut = Stokes::PFtotal;
	exprName = String("ComplexFractionalTotalPolarization");
      } else if (tA==Stokes::PFlinear) {
	typeOut = Stokes::PFlinear;
	exprName = String("ComplexFractionalLinearPolarization");
      } else {
	*itsLog << "Cannot form Complex image for this amplitude image"
		<< endl;
	*itsLog << "Expect linear, total, or fractional polarization"
		<< LogIO::EXCEPTION;
      }
      //
      LatticeExprNode node0(2.0*fac*pp);
      LatticeExprNode node(formComplex(aa*cos(node0),aa*sin(node0)));
      LatticeExpr<Complex> le(node);
      pExpr = new ImageExpr<Complex>(le, exprName);
      fiddleStokesCoordinate(*pExpr, typeOut);
    }

    // Copy to output
    pOutComplex->setCoordinateInfo(pExpr->coordinates());
    LatticeUtilities::copyDataAndMask(*itsLog, *pOutComplex, *pExpr);
    //
    auto p = itsImPol->imageInterface();
    copyMiscellaneous (*pOutComplex, *p);
    //
    delete pExpr; pExpr = 0;
    delete pOutComplex; pOutComplex = 0;
  }


  // Private functions

  Bool ImagePol::copyImage(ImageInterface<Float>*& out, 
			   const ImageInterface<Float>& in, 
			   const String& outfile, Bool overwrite){

    // if no outfile, just make out image from the input image
    if (outfile.empty()){
      out = in.cloneII();
      return true;
    }

    // The user wants to write the image out; verify file
    if (!overwrite) {
       NewFile validfile;
       String errmsg;
       if (!validfile.valueOK(outfile, errmsg)) {
           *itsLog << errmsg << LogIO::EXCEPTION;
       }
    }

    // Create output image
    out= new PagedImage<Float>(in.shape(), in.coordinates(), outfile);
    if (out == 0) {
      *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
    } else {
      *itsLog << LogIO::NORMAL << "Creating image '" << outfile
	      << "' of shape " << out->shape() << LogIO::POST;
    }

    // Make mask
    if (in.isMasked()) makeMask(*out, false);

    // Copy stuff
    LatticeUtilities::copyDataAndMask (*itsLog, *out, in);
    ImageUtilities::copyMiscellaneous (*out, in);
    
    return true;

  }

  Bool ImagePol::makeMask(ImageInterface<Float>& out, Bool init){

    if (out.canDefineRegion()) {

      // Generate mask name if not given
      String maskName = out.makeUniqueRegionName(String("mask"), 0);
      
      // Make the mask if it does not exist
      if (!out.hasRegion (maskName, RegionHandler::Masks)) {
         out.makeMask (maskName, true, true, init, true);
	 /* if (init) {
            *itsLog << LogIO::NORMAL << "Created and initialized mask `" << maskName << "'" << LogIO::POST;
         } else {
            *itsLog << LogIO::NORMAL << "Created mask `" << maskName << "'" << LogIO::POST;
	    }*/
      }
      return true;
    } else {
      *itsLog << LogIO::WARN
	      << "Cannot make requested mask for this type of image" << endl;
      return false;
    }
    
  }
  Bool ImagePol::makeMask(ImageInterface<Complex>& out, Bool init){

    if (out.canDefineRegion()) {

      // Generate mask name if not given
      String maskName = out.makeUniqueRegionName(String("mask"), 0);
      
      // Make the mask if it does not exist
      if (!out.hasRegion (maskName, RegionHandler::Masks)) {
	out.makeMask (maskName, true, true, init, true);
	/* if (init) {
	 *itsLog << LogIO::NORMAL << "Created and initialized mask `" << maskName << "'" << LogIO::POST;
         } else {
	 *itsLog << LogIO::NORMAL << "Created mask `" << maskName << "'" << LogIO::POST;
	 }*/
      }
      return true;
    } else {
      *itsLog << LogIO::WARN
	      << "Cannot make requested mask for this type of image" << endl;
      return false;
    }
    
  }


  Bool ImagePol::makeImage(ImageInterface<Float>*& out, 
			   const String& outfile, const CoordinateSystem& cSys,
			   const IPosition& shape, Bool isMasked,
			   Bool tempAllowed) {
    // Verify outfile
    if (outfile.empty()) {
      if (!tempAllowed) return false;
    }
    // else {
    //  NewFile validfile;
    //  String errmsg;
    //  if (!validfile.valueOK(outfile, errmsg)) {
    //	*itsLog << errmsg << LogIO::EXCEPTION;
    //  }
    //}


    uInt ndim = shape.nelements();
    if (ndim != cSys.nPixelAxes()) {
      *itsLog << "Supplied CoordinateSystem and image shape are inconsistent"
	      << LogIO::EXCEPTION;
    }
    if (outfile.empty()) {
       out = new TempImage<Float>(shape, cSys);
       if (out == 0) {
          *itsLog << "Failed to create TempImage" << LogIO::EXCEPTION;
       }
       *itsLog << LogIO::NORMAL << "Creating (temp)image of shape "
	       << out->shape() << LogIO::POST;
    } else {
       out = new PagedImage<Float>(shape, cSys, outfile);
       if (out == 0) {
	 *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
       }
       *itsLog << LogIO::NORMAL << "Creating image '" << outfile
	       << "' of shape " << out->shape() << LogIO::POST;
    }
    //
    if (isMasked) {
       makeMask(*out, true);
    }

    return true;
  }

Bool ImagePol::_makeImage(
    ImageInterface<Complex>*& out, const String& outfile,
    const CoordinateSystem& cSys, const IPosition& shape, Bool isMasked,
    Bool tempAllowed
) {
    const auto ndim = shape.nelements();
    if (ndim != cSys.nPixelAxes()) {
        *itsLog << "Supplied CoordinateSystem and image shape are inconsistent"
	        << LogIO::EXCEPTION;
    }
    if (outfile.empty()) {
        if (tempAllowed) {
            out = new TempImage<Complex>(shape, cSys);
            if (! out) {
                *itsLog << "Failed to create TempImage" << LogIO::EXCEPTION;
            }
            *itsLog << LogIO::NORMAL << "Creating (temp)image of shape "
	            << out->shape() << LogIO::POST;
        }
        else {
            return false;
        }
    }
    else {
        out = new PagedImage<Complex>(shape, cSys, outfile);
        if (! out) {
	        *itsLog << "Failed to create PagedImage" << LogIO::EXCEPTION;
        }
        *itsLog << LogIO::NORMAL << "Creating image '" << outfile
	        << "' of shape " << out->shape() << LogIO::POST;
    }
    if (isMasked) {
        makeMask(*out, true);
    }
    return true;
}

SPIIC ImagePol::_makeImage( 
    const String& outfile, const CoordinateSystem& cSys,
    const IPosition& shape, bool isMasked,
    bool tempAllowed
) {
    ImageInterface<Complex>* out = NULL;
    const auto res = _makeImage(
        out, outfile, cSys, shape, isMasked, tempAllowed
    );
    if (res) {
        return shared_ptr<ImageInterface<Complex>>(out);
    }
    else {
        if (out) {
            delete out;
        }
        return shared_ptr<ImageInterface<Complex>>();
    }
}

  void ImagePol::copyMiscellaneous (ImageInterface<Complex>& out,
				    const ImageInterface<Float>& in)
  {
    out.setMiscInfo(in.miscInfo());
    out.appendLog(in.logger());
  }

  void ImagePol::copyMiscellaneous (ImageInterface<Float>& out,
				    const ImageInterface<Float>& in)
  {
    out.setMiscInfo(in.miscInfo());
    out.appendLog(in.logger());
  }


  void ImagePol::fiddleStokesCoordinate(ImageInterface<Float>& ie,
					Stokes::StokesTypes type)
  {
    CoordinateSystem cSys = ie.coordinates();
    //
    Int afterCoord = -1;
    Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
    //
    Vector<Int> which(1);
    which(0) = Int(type);
    StokesCoordinate stokes(which);
    cSys.replaceCoordinate(stokes, iStokes);
    ie.setCoordinateInfo(cSys);
  }

  void ImagePol::fiddleStokesCoordinate(ImageInterface<Complex>& ie, 
					Stokes::StokesTypes type)
  {
    CoordinateSystem cSys = ie.coordinates();
    //
    Int afterCoord = -1;
    Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
    //
    Vector<Int> which(1);
    which(0) = Int(type);
    StokesCoordinate stokes(which);
    cSys.replaceCoordinate(stokes, iStokes);
    ie.setCoordinateInfo(cSys);
  }


  Bool ImagePol::makeIQUVImage (ImageInterface<Float>*& pImOut, 
				const String& outfile, 
				Double sigma, Double pa0, 
				const Vector<Float>& rm, 
				const IPosition& shape,
				Double f0, Double dF)
    //
    // Must be 4D
    //
  {
    AlwaysAssert(shape.nelements()==4,AipsError);
    AlwaysAssert(shape(2)==4,AipsError);
    //
    CoordinateSystem cSys;
    CoordinateUtil::addDirAxes(cSys);
    //
    Vector<Int> whichStokes(4);
    whichStokes(0) = Stokes::I;
    whichStokes(1) = Stokes::Q;
    whichStokes(2) = Stokes::U;
    whichStokes(3) = Stokes::V;
    StokesCoordinate stokesCoord(whichStokes);
    cSys.addCoordinate(stokesCoord);
    //
    const Int nchan = shape(3);
    Double df = dF / nchan;
    Double refpix = 0.0;
    SpectralCoordinate spectCoord(MFrequency::TOPO, f0, df, refpix, f0);
    cSys.addCoordinate(spectCoord);
   
    // Centre reference pixel
    centreRefPix (cSys, shape);
   
    // Make image 
    makeImage (pImOut, outfile, cSys, shape, false, true);
    //
    uInt stokesAxis = 2;
    uInt spectralAxis = 3;
   
    // Fill image with I, Q, U and V. 
    fillIQUV (*pImOut, stokesAxis, spectralAxis, rm, pa0);

    // Add noise 
    Array<Float> slice = pImOut->get();
    Float maxVal = max(slice);
    Float t = sigma * maxVal;
    *itsLog << LogIO::NORMAL << "Using sigma            = "
	    << t << LogIO::POST;
    MLCG gen;
    Normal noiseGen(&gen, 0.0, t*t);
    addNoise(slice, noiseGen);
    pImOut->put(slice);
    return true;
  }

Bool  ImagePol::fillIQUV (ImageInterface<Float>& im, uInt stokesAxis, 
                         uInt spectralAxis, const Vector<Float>& rm, 
                         Float pa0)
  //
  // Image must be 4D
  //
{

  // Find spectral coordinate
  const CoordinateSystem& cSys = im.coordinates();   
  Int spectralCoord, iDum;
  cSys.findPixelAxis(spectralCoord, iDum, spectralAxis);
  const SpectralCoordinate& sC = cSys.spectralCoordinate(spectralCoord);
  //
  IPosition shape = im.shape();
  Double c = QC::c( ).getValue(Unit("m/s"));
  Double lambdasq;
  MFrequency freq;
  IPosition blc(4,0);
  IPosition trc(shape-1);
  //
  Double ii = 2.0;                      // arbitrary
  Double vv = 0.05 * ii;                // arbitrary
  const Int n = shape(3);
  const uInt nRM = rm.nelements();
  for (Int i=0; i<n; i++) {
    if (!sC.toWorld(freq, Double(i))) {
      *itsLog << sC.errorMessage() << LogIO::EXCEPTION;
    }
    Double fac = c / freq.get(Unit("Hz")).getValue();
    lambdasq = fac*fac;
    //
    Double chi = rm(0)*lambdasq + pa0;
    Double q = cos(2*chi);
    Double u = sin(2*chi);
    //
    if (nRM > 1) {
      for (uInt j=1; j<nRM; j++) {
	chi = rm(j)*lambdasq + pa0;
	q += cos(2*chi);
	u += sin(2*chi);
      }
    }
    q = q / Double(nRM);
    u = u / Double(nRM);
    //
    blc(spectralAxis) = i;                // channel
    trc(spectralAxis) = i;
    {
      blc(stokesAxis) = 1;                // Q       
      trc(stokesAxis) = 1;                
      Slicer sl(blc, trc, Slicer::endIsLast);
      SubImage<Float> subImage(im, sl, true);
      subImage.set(q);
    }
    {
      blc(stokesAxis) = 2;                // U
      trc(stokesAxis) = 2;                
      Slicer sl(blc, trc, Slicer::endIsLast);
      SubImage<Float> subImage(im, sl, true);
      subImage.set(u);
    }
  }
  // 
  blc(spectralAxis) = 0;  
  trc(spectralAxis) = n-1;
  {
    blc(stokesAxis) = 0;                // I
    trc(stokesAxis) = 0;                
    Slicer sl(blc, trc, Slicer::endIsLast);
    SubImage<Float> subImage(im, sl, true);
    subImage.set(ii);
  }
  {
    blc(stokesAxis) = 3;                // V
    trc(stokesAxis) = 3;                
    Slicer sl(blc, trc, Slicer::endIsLast);
    SubImage<Float> subImage(im, sl, true);
    subImage.set(vv);
  }
  return true;
  
}

void ImagePol::addNoise (Array<Float>& slice, Normal& noiseGen) 
{
   Bool deleteIt;
   Float* p = slice.getStorage(deleteIt);
   for (uInt i=0; i<slice.nelements(); i++) {
      p[i] += noiseGen();
   }
   slice.putStorage(p, deleteIt);
}

void ImagePol::centreRefPix (CoordinateSystem& cSys, const IPosition& shape)
{
   Int after = -1;
   Int iS = cSys.findCoordinate(Coordinate::STOKES, after);
   Int sP = -1;
   if (iS>=0) {
      Vector<Int> pixelAxes = cSys.pixelAxes(iS);
      sP = pixelAxes(0);
   }
   Vector<Double> refPix = cSys.referencePixel();
   for (Int i=0; i<Int(refPix.nelements()); i++) {
     if (i!=sP) refPix(i) = Double(shape(i) / 2);
   }     
   cSys.setReferencePixel(refPix);
}


Stokes::StokesTypes ImagePol::stokesType(const CoordinateSystem& cSys) 
{
  Stokes::StokesTypes type = Stokes::Undefined;
  Int afterCoord = -1;
  Int iStokes = cSys.findCoordinate(Coordinate::STOKES, afterCoord);
  if (iStokes >=0) {
    Vector<Int> which = cSys.stokesCoordinate(iStokes).stokes();
    if (which.nelements()>1) {
        *itsLog << "Stokes axis must be of length unity" << LogIO::EXCEPTION;
      } else {
         type = Stokes::type(which(0));
      }
  } else {
    *itsLog << "No StokesCoordinate" << LogIO::EXCEPTION;
   }
   return type;
}


} // end of  casa namespace