lyotCoronagraph.hpp

Go to the documentation of this file.

/** \file lyotCoronagraph.hpp
 * \brief Declares and defines a class to describe and optimize a Lyot Coronagraph.
 * \ingroup imaging_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017, 2018 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef __lyotCoronagraph_hpp__
#define __lyotCoronagraph_hpp__
 
#include <iostream>
 
#include "imagingArray.hpp"
#include "imagingUtils.hpp"
#include "fraunhoferPropagator.hpp"
 
#include "../ioutils/fits/fitsFile.hpp"
#include "../ioutils/fits/fitsUtils.hpp"
 
#include "../error/mxErrorOld.hpp"
 
#include "../improc/eigenImage.hpp"
#include "../improc/eigenCube.hpp"
 
namespace mx
{
 
namespace wfp
{
 
/// The Lyot Coronagraph
/** A generalized Lyot coronagraph, which can include a pupil apodization, a focal plane mask with complex transmission,
 * and a Lyot stop.  Light is propagated through the coronagraph plane-by-plane, and the complex wavefront can be
 * accessed at any plane. Also provides functions for optimizing the pupil apodizer and focal plane mask.
 *
 * \ingroup coronagraphs
 */
template <typename _realT, typename _fpmaskFloatT>
struct lyotCoronagraph
{
  public:
    /// The real floating point type
    typedef _realT realT;
 
    /// The real floating point type for mask calculations.
    typedef _fpmaskFloatT fpmaskFloatT;
 
    /// The complex floating point type
    typedef std::complex<realT> complexT;
 
    /// The wavefront complex field type
    typedef mx::wfp::imagingArray<std::complex<realT>, fftwAllocator<std::complex<realT>>, 0> complexFieldT;
 
    /// The focal plane mask type
    typedef Eigen::Array<std::complex<fpmaskFloatT>, Eigen::Dynamic, Eigen::Dynamic> fpMaskT;
 
    /// The image type
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
 
  public:
    /// The directory where coronagraph files are stored
    std::string m_fileDir{ "coron" };
 
    /// The linear size of the wavefront in pixels
    int m_wfSz{ 0 };
 
    /// Image containing the pupil apodization.
    imageT m_pupilApodizer;
 
    int m_maskSource;       ///< 0= read from file, 1 = constructed by makeFocalMask, 2 = trans optimized.
    std::string m_maskFile; ///< Name of file from which mask was loaded.
 
    realT m_maskRad;   ///< Radius of mask if it was constructed.
    realT m_maskTrans; ///< Transmission of mask if it was constructed.
 
    /// The focal plane  mask.
    fpMaskT m_focalMask;
 
    /// Image containing the lyot stop.
    imageT m_lyotStop;
 
    complexFieldT m_focalPlane;
 
    bool m_savePreMaskFocalPlane{ false };
    complexFieldT m_preMaskFocalPlane;
 
    bool m_savePreLyotPupilPlane{ false };
    complexFieldT m_preLyotPupilPlane;
 
    /// Fraunhofer propagator
    fraunhoferPropagator<complexFieldT> m_fp;
 
  public:
    /// Default c'tor
    lyotCoronagraph();
 
    /// Get the wavefront size in pixels
    /**
     * \returns the wavefront size in pixels
     */
    int wfSz();
 
    /// Set the wavefront size in pixels.
    /**
     */
    void wfSz( int sz /**< [in] is the new size */ );
 
    /// Make the focal plane mask
    void makeFocalMask( realT rad, fpmaskFloatT trans = 0.0, int sz = 0.0 );
 
    /// Load the apodizer from a FITS file
    /**
     * \returns 0 on success.
     * \returns -1 on error.
     */
    int loadApodizer( const std::string &apodName /**< [in] is the name of the FITS file containing the apodizer. */ );
 
    /// Load the focal plane mask from a FITS file
    /**
     * \returns 0 on success.
     * \returns -1 on error.
     */
    int loadFocalMask(
        const std::string &fpmName /**< [in] is the name of the FITS file containing the focal plane mask. */ );
 
    /// Load the Lyot stop from a FITS file
    /**
     * \returns 0 on success.
     * \returns -1 on error.
     */
    int loadLyotStop( const std::string &lyotName /**< [in] is the name of the FITS file containing the lyot stop. */ );
 
    /// Load the components of the coronagraph from FITS files
    /**
     * \returns 0 on success.
     * \returns -1 on error.
     */
    int
    loadCoronagraph( const std::string &apodName, ///< [in] is the name of the FITS file containing the apodizer.
                     const std::string &fpmName, ///< [in] is the name of the FITS file containing the focal plane mask.
                     const std::string &lyotName ///< [in] is the name of the FITS file containing the lyot stop.
    );
 
    /// Load the components of the coronagraph based in its base name
    /** Looks in _filDir for the files.
     *
     * \returns 0 on success.
     * \returns -1 on error.
     */
    int
    loadCoronagraph( const std::string
                         &cName /**< [in] is the base name of the coronagraph without directory or file extension. */ );
 
    void applyApodizer( complexFieldT &pupilPlane );
    void applyFocalMask( complexFieldT &focalPlane );
    void applyLyotStop( complexFieldT &lyotPlane );
 
    /// Propagate the given pupil-plane wavefront through the coronagraph.
    int propagate(
        complexFieldT
            &pupilPlane /**< [in.out] The wavefront at the input pupil plane.  It is modified by the coronagraph. */ );
 
    int
    propagate( imageT &fpIntensity, ///< [out] The intensity image in the focal plane.  This should be pre-allocated.
               complexFieldT
                   &pupilPlane ///< [in.out] The wavefront at the input pupil plane.  It is modified by the coronagraph.
    );
 
    /// Propagate the given pupil-plane wavefront without the coronagraph.
    /** For a Lyot coronagraph, this applies the pupil apodization and Lyot stop, but not the FPM, to the given
     * pupil-plane wavefront such that the result will produce the non-coronagraphic (off-axis) PSF.
     */
    int propagateNC(
        complexFieldT
            &pupilPlane /**< [in.out] The wavefront at the input pupil plane.  It is modified by the coronagraph. */ );
 
    int propagateNC(
        imageT &fpIntensity, ///< [out] The intensity image in the focal plane.  This should be pre-allocated.
        complexFieldT
            &pupilPlane ///< [in.out] The wavefront at the input pupil plane.  It is modified by the coronagraph.
    );
 
    /// Optimize the pupil amplitude apodization and focal-plane mask complex transmission.
    /** Uses the algorithm in Guyon (2014) \cite guyon_2014 for optimizing an Apodized-Pupil Lyot Complex Mask
     * Coronagraph (APLCMC).
     *
     */
    void
    optimizeAPLCMC( imageT &geomPupil,       ///< The geometric pupil mask, binary 1/0 transmission.
                    realT fpmRadPix,         ///< The radius in pixels of the FPM.
                    realT relTol,            ///< Relative tolerance for convergence.
                    realT absTol,            ///< Absolute tolerance for convergence.
                    int maxIter,             ///< Maximum number of iterations to allow.
                    const std::string &cname ///< Name of coronagraph, used as base-name for output files in m_fileDir.
    );
 
    void
    optimizeAPLCMC( imageT &geomPupil, ///< The geometric pupil mask, binary 1/0 transmission.
                    realT fpmRadPix,   ///< The radius in pixels of the FPM.
                    imageT &fpmPhase,
                    realT relTol,            ///< Relative tolerance for convergence.
                    realT absTol,            ///< Absolute tolerance for convergence.
                    int maxIter,             ///< Maximum number of iterations to allow.
                    const std::string &cname ///< Name of coronagraph, used as base-name for output files in m_fileDir.
    );
};
 
template <typename _realT, typename _fpmaskFloatT>
lyotCoronagraph<_realT, _fpmaskFloatT>::lyotCoronagraph()
{
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::wfSz()
{
    return m_wfSz;
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::wfSz( int sz )
{
    m_wfSz = sz;
 
    m_fp.setWavefrontSizePixels( sz );
    m_focalPlane.resize( sz, sz );
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::makeFocalMask( _realT rad, _fpmaskFloatT trans, int sz )
{
    if( sz == 0 )
    {
        sz = 2 * rad + 1;
        if( sz % 2 == 1 )
            ++sz;
    }
 
    // Make a real circular mask
    improc::eigenImage<fpmaskFloatT> fpm( sz, sz );
    mx::wfp::circularPupil( fpm, 0., rad );
 
    // Convert to complex amplitude
    m_focalMask.resize( sz, sz );
    for( int i = 0; i < sz; ++i )
    {
        for( int j = 0; j < sz; ++j )
        {
            if( fpm( i, j ) == 1 )
                m_focalMask( i, j ) = std::complex<fpmaskFloatT>( trans, 0 );
            else
                m_focalMask( i, j ) = std::complex<fpmaskFloatT>( 1, 0 );
        }
    }
 
    m_maskSource = 1;
    m_maskFile = "";
    m_maskRad = rad;
    m_maskTrans = 0.0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::loadApodizer( const std::string &apodName )
{
    fits::fitsFile<_realT> ff;
 
    if( !!ff.read( m_pupilApodizer, apodName ) )
        return -1;
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::loadFocalMask( const std::string &fpmName )
{
    fits::fitsFile<_realT> ff;
    improc::eigenCube<fpmaskFloatT> fpm;
 
    if( !!ff.read( fpm, fpmName ) )
        return -1;
 
    if( fpm.planes() == 1 )
    {
        m_focalMask = fpm.image( 0 );
    }
    else if( fpm.planes() == 2 )
    {
        m_focalMask.resize( fpm.rows(), fpm.cols() );
 
        for( int r = 0; r < fpm.rows(); ++r )
        {
            for( int c = 0; c < fpm.cols(); ++c )
            {
                m_focalMask( r, c ) =
                    fpm.image( 0 )( r, c ) * exp( std::complex<fpmaskFloatT>( 0, fpm.image( 1 )( r, c ) ) );
            }
        }
    }
    else
    {
        std::cerr << "too many planes in focal mask file\n";
        return -1;
    }
 
    m_maskSource = 0;
    m_maskFile = fpmName;
    m_maskRad = 0.0;
    m_maskTrans = 0.0;
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::loadLyotStop( const std::string &lyotName )
{
    fits::fitsFile<_realT> ff;
 
    if( !!ff.read( m_lyotStop, lyotName ) )
        return -1;
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::loadCoronagraph( const std::string &apodName,
                                                             const std::string &fpmName,
                                                             const std::string &lyotName )
{
    if( loadApodizer( apodName ) < 0 )
        return -1;
    if( loadFocalMask( fpmName ) < 0 )
        return -1;
    if( loadLyotStop( lyotName ) < 0 )
        return -1;
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::loadCoronagraph( const std::string &cName )
{
    if( m_fileDir == "" )
    {
        mxError( "lyotCoronagraph", MXE_PARAMNOTSET, "file directory (fileDir) not set." );
        return -1;
    }
 
    std::string apodName = m_fileDir + "/" + cName + "_apod.fits";
    std::string fpmName = m_fileDir + "/" + cName + "_fpm.fits";
    std::string lyotName = m_fileDir + "/" + cName + "_lyot.fits";
 
    return loadCoronagraph( apodName, fpmName, lyotName );
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::applyApodizer( complexFieldT &pupilPlane )
{
    int sz = m_pupilApodizer.rows();
 
    realT w = 0.5 * ( sz - 1.0 );
 
    realT xc = 0.5 * ( pupilPlane.rows() - 1 );
    realT yc = 0.5 * ( pupilPlane.cols() - 1 );
 
    for( int i = 0; i < pupilPlane.rows(); ++i )
    {
        for( int j = 0; j < pupilPlane.cols(); ++j )
        {
            if( i >= xc - w && i <= xc + w && j >= yc - w && j <= yc + w )
                pupilPlane( i, j ) *= m_pupilApodizer( (int)( i - ( xc - w ) ), (int)( j - ( yc - w ) ) );
            else
                pupilPlane( i, j ) *= 0;
        }
    }
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::applyFocalMask( complexFieldT &focalPlane )
{
    int sz = m_focalMask.rows();
 
    realT w = 0.5 * ( sz - 1.0 );
 
    realT xc = 0.5 * ( focalPlane.rows() - 1 );
    realT yc = 0.5 * ( focalPlane.cols() - 1 );
 
    for( int i = 0; i < sz; ++i )
    {
        for( int j = 0; j < sz; ++j )
        {
            focalPlane( xc - w + i, yc - w + j ) *= m_focalMask( i, j );
        }
    }
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::applyLyotStop( complexFieldT &lyotPlane )
{
    int sz = m_lyotStop.rows();
 
    realT w = 0.5 * ( sz - 1.0 );
 
    realT xc = 0.5 * ( lyotPlane.rows() - 1 );
    realT yc = 0.5 * ( lyotPlane.cols() - 1 );
 
    for( int i = 0; i < lyotPlane.rows(); ++i )
    {
        for( int j = 0; j < lyotPlane.cols(); ++j )
        {
            if( i >= xc - w && i <= xc + w && j >= yc - w && j <= yc + w )
                lyotPlane( i, j ) *= m_lyotStop( (int)( i - ( xc - w ) ), (int)( j - ( yc - w ) ) );
            else
                lyotPlane( i, j ) = 0;
        }
    }
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::propagate( complexFieldT &pupilPlane )
{
    applyApodizer( pupilPlane );
 
    m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
 
    if( m_savePreMaskFocalPlane )
    {
        m_preMaskFocalPlane = m_focalPlane;
    }
 
    applyFocalMask( m_focalPlane );
 
    m_fp.propagateFocalToPupil( pupilPlane, m_focalPlane );
 
    if( m_savePreLyotPupilPlane )
    {
        m_preLyotPupilPlane = pupilPlane;
    }
 
    applyLyotStop( pupilPlane );
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::propagate( imageT &fpIntensity, complexFieldT &pupilPlane )
{
    propagate( pupilPlane );
 
    m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
 
    int x0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
    int y0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
 
    extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
 
    return 0;
}
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::propagateNC( complexFieldT &pupilPlane )
{
    applyApodizer( pupilPlane );
 
    m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
 
    m_fp.propagateFocalToPupil( pupilPlane, m_focalPlane );
 
    applyLyotStop( pupilPlane );
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
int lyotCoronagraph<_realT, _fpmaskFloatT>::propagateNC( imageT &fpIntensity, complexFieldT &pupilPlane )
{
    propagateNC( pupilPlane );
 
    m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
 
    int x0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
    int y0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
 
    extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
 
    return 0;
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::optimizeAPLCMC(
    imageT &geomPupil, realT fpmRadPix, realT relTol, realT absTol, int maxIter, const std::string &cname )
{
    complexFieldT focalPlane, pupilPlane;
 
    pupilPlane.resize( m_wfSz, m_wfSz );
    focalPlane.resize( m_wfSz, m_wfSz );
 
    mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> mask( m_wfSz, m_wfSz );
    mx::wfp::circularPupil( mask, 0., fpmRadPix );
    /*for(int c=0; c< mask.cols()*0.5; ++c)
    {
       for(int r=0; r < mask.rows(); ++r) mask(r,c) = 0;
    }*/
 
    // Initialize pupilImage
    mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> pupilImage( m_wfSz, m_wfSz );
    pupilImage.setZero();
 
    int gpLLi = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.rows() - 1 );
    int gpLLj = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.cols() - 1 );
 
    int gpURi = gpLLi + geomPupil.rows();
    int gpURj = gpLLj + geomPupil.cols();
 
    for( int i = 0; i < geomPupil.rows(); ++i )
    {
        for( int j = 0; j < geomPupil.cols(); ++j )
        {
            pupilImage( gpLLi + i, gpLLj + j ) = geomPupil( i, j );
        }
    }
 
    realT lastLambdaA, LambdaA;
 
    lastLambdaA = 1;
    int n;
    std::string reason;
    for( n = 0; n < maxIter; ++n )
    {
        mx::wfp::makeComplexPupil( pupilPlane, pupilImage, m_wfSz );
 
        m_fp.propagatePupilToFocal( focalPlane, pupilPlane );
 
        for( int i = 0; i < m_wfSz; ++i )
        {
            for( int j = 0; j < m_wfSz; ++j )
            {
                focalPlane( i, j ) *= mask( i, j );
            }
        }
 
        m_fp.propagateFocalToPupil( pupilPlane, focalPlane );
 
        for( int i = 0; i < m_wfSz; ++i )
            for( int j = 0; j < m_wfSz; ++j )
                pupilImage( i, j ) = abs( pupilPlane( i, j ) );
 
        LambdaA = 0;
        for( int i = 0; i < m_wfSz; ++i )
        {
            for( int j = 0; j < m_wfSz; ++j )
            {
                if( i >= gpLLi && i < gpURi && j >= gpLLj && j < gpURj )
                {
                    pupilImage( i, j ) *= geomPupil( i - gpLLi, j - gpLLj );
                    if( pupilImage( i, j ) > LambdaA )
                        LambdaA = pupilImage( i, j );
                }
                else
                {
                    pupilImage( i, j ) = 0;
                }
            }
        }
        // LambdaA = 1.0/LambdaA;
        for( int i = 0; i < m_wfSz; ++i )
            for( int j = 0; j < m_wfSz; ++j )
                pupilImage( i, j ) /= LambdaA;
 
        std::cout << n << " " << LambdaA << "\n";
        if( fabs( LambdaA - lastLambdaA ) < absTol )
        {
            std::cout << "Converged on absTol.\n";
            reason = "absTol";
            break;
        }
        if( fabs( ( LambdaA - lastLambdaA ) / lastLambdaA ) < relTol )
        {
            std::cout << "Converged on relTol.\n";
            reason = "relTol";
            break;
        }
 
        if( n == maxIter - 1 )
        {
            std::cout << "maxIter reached.\n";
            reason = "maxIter";
        }
 
        lastLambdaA = LambdaA;
    }
    if( LambdaA > 1 )
        LambdaA = 1;
 
    std::cout << "LambdaA: = " << LambdaA << "\n";
 
    realT trans = 1.0 - 1.0 / LambdaA;
 
    int pupSize = geomPupil.rows();
 
    m_pupilApodizer.resize( pupSize, pupSize );
 
    extractBlock( m_pupilApodizer,
                  0,
                  pupSize,
                  0,
                  pupSize,
                  pupilImage,
                  0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ),
                  0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ) );
 
    makeFocalMask( fpmRadPix, trans, pupSize );
 
    /*for(int c=0; c< m_focalMask.cols()*0.5; ++c)
    {
       for(int r=0; r < m_focalMask.rows(); ++r)
       {
          if(m_focalMask(r,c).real()==1) m_focalMask(r,c) = 0;
       }
    }*/
 
    m_maskSource = 2;
 
    m_lyotStop = geomPupil;
 
    fits::fitsHeader head;
 
    head.append( "", fits::fitsCommentType(), "----------------------------------------" );
    head.append( "", fits::fitsCommentType(), "lyotCoronagraph optimization Parameters:" );
    head.append( "", fits::fitsCommentType(), "----------------------------------------" );
    head.append<int>( "WFSZ", m_wfSz, "Size of wavefront used for FFTs (pixels)" );
    head.append<realT>( "FPMRADPX", fpmRadPix, "input radius of focal plane mask (pixels)" );
    head.append<realT>( "ABSTOL", absTol, "input absolute tolerance" );
    head.append<realT>( "RELTOL", relTol, "input relative tolerance" );
    head.append<int>( "MAXITER", maxIter, "input maximum iterations" );
    head.append<int>( "NITER", n, "actual number of iterations" );
    head.append<std::string>( "XREASON", reason, "reason for convergence" );
    head.append<realT>( "FPMTRANS", trans, "transmission of FPM" );
 
    fits::fitsFile<double> ff;
 
    std::string fname = "coron/" + cname + "_apod.fits";
 
    ff.write( fname, m_pupilApodizer, head );
 
    fname = "coron/" + cname + "_fpm.fits";
    improc::eigenImage<fpmaskFloatT> fpm( m_focalMask.rows(), m_focalMask.cols() );
    for( int r = 0; r < m_focalMask.rows(); ++r )
    {
        for( int c = 0; c < m_focalMask.cols(); ++c )
        {
            fpm( r, c ) = m_focalMask( r, c ).real();
        }
    }
 
    ff.write( fname, fpm, head );
 
    fname = "coron/" + cname + "_lyot.fits";
    ff.write( fname, m_lyotStop, head );
}
 
template <typename _realT, typename _fpmaskFloatT>
void lyotCoronagraph<_realT, _fpmaskFloatT>::optimizeAPLCMC( imageT &geomPupil,
                                                             realT fpmRadPix,
                                                             imageT &fpmPhase,
                                                             realT relTol,
                                                             realT absTol,
                                                             int maxIter,
                                                             const std::string &cname )
{
    complexFieldT focalPlane, pupilPlane;
 
    pupilPlane.resize( m_wfSz, m_wfSz );
    focalPlane.resize( m_wfSz, m_wfSz );
 
    mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> mask( m_wfSz, m_wfSz );
    mx::wfp::circularPupil( mask, 0., fpmRadPix );
 
    // Initialize pupilImage
    mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> pupilImage( m_wfSz, m_wfSz );
    pupilImage.setZero();
 
    int gpLLi = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.rows() - 1 );
    int gpLLj = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.cols() - 1 );
 
    int gpURi = gpLLi + geomPupil.rows();
    int gpURj = gpLLj + geomPupil.cols();
 
    for( int i = 0; i < geomPupil.rows(); ++i )
    {
        for( int j = 0; j < geomPupil.cols(); ++j )
        {
            pupilImage( gpLLi + i, gpLLj + j ) = geomPupil( i, j );
        }
    }
 
    realT lastLambdaA, LambdaA;
 
    lastLambdaA = 1;
    int n;
    std::string reason;
    for( n = 0; n < maxIter; ++n )
    {
        mx::wfp::makeComplexPupil( pupilPlane, pupilImage, m_wfSz );
 
        m_fp.propagatePupilToFocal( focalPlane, pupilPlane );
 
        for( int i = 0; i < m_wfSz; ++i )
        {
            for( int j = 0; j < m_wfSz; ++j )
            {
                focalPlane( i, j ) *= mask( i, j ) * exp( std::complex<realT>( 0, fpmPhase( i, j ) ) );
                ;
            }
        }
 
        m_fp.propagateFocalToPupil( pupilPlane, focalPlane );
 
        for( int i = 0; i < m_wfSz; ++i )
            for( int j = 0; j < m_wfSz; ++j )
                pupilImage( i, j ) = abs( pupilPlane( i, j ) );
 
        LambdaA = 0;
        for( int i = 0; i < m_wfSz; ++i )
        {
            for( int j = 0; j < m_wfSz; ++j )
            {
                if( i >= gpLLi && i < gpURi && j >= gpLLj && j < gpURj )
                {
                    pupilImage( i, j ) *= geomPupil( i - gpLLi, j - gpLLj );
                    if( pupilImage( i, j ) > LambdaA )
                        LambdaA = pupilImage( i, j );
                }
                else
                {
                    pupilImage( i, j ) = 0;
                }
            }
        }
        // LambdaA = 1.0/LambdaA;
        for( int i = 0; i < m_wfSz; ++i )
            for( int j = 0; j < m_wfSz; ++j )
                pupilImage( i, j ) /= LambdaA;
 
        std::cout << n << " " << LambdaA << "\n";
        if( fabs( LambdaA - lastLambdaA ) < absTol )
        {
            std::cout << "Converged on absTol.\n";
            reason = "absTol";
            break;
        }
        if( fabs( ( LambdaA - lastLambdaA ) / lastLambdaA ) < relTol )
        {
            std::cout << "Converged on relTol.\n";
            reason = "relTol";
            break;
        }
 
        if( n == maxIter - 1 )
        {
            std::cout << "maxIter reached.\n";
            reason = "maxIter";
        }
 
        lastLambdaA = LambdaA;
    }
    if( LambdaA > 1 )
        LambdaA = 1;
 
    std::cout << "LambdaA: = " << LambdaA << "\n";
 
    realT trans = 1.0 - 1.0 / LambdaA;
 
    int pupSize = geomPupil.rows();
 
    m_pupilApodizer.resize( pupSize, pupSize );
 
    extractBlock( m_pupilApodizer,
                  0,
                  pupSize,
                  0,
                  pupSize,
                  pupilImage,
                  0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ),
                  0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ) );
 
    makeFocalMask( fpmRadPix, trans, pupSize );
    m_maskSource = 2;
 
    m_lyotStop = geomPupil;
 
    fits::fitsHeader head;
 
    head.append( "", fits::fitsCommentType(), "----------------------------------------" );
    head.append( "", fits::fitsCommentType(), "lyotCoronagraph optimization Parameters:" );
    head.append( "", fits::fitsCommentType(), "----------------------------------------" );
    head.append<int>( "WFSZ", m_wfSz, "Size of wavefront used for FFTs (pixels)" );
    head.append<realT>( "FPMRADPX", fpmRadPix, "input radius of focal plane mask (pixels)" );
    head.append<realT>( "ABSTOL", absTol, "input absolute tolerance" );
    head.append<realT>( "RELTOL", relTol, "input relative tolerance" );
    head.append<int>( "MAXITER", maxIter, "input maximum iterations" );
    head.append<int>( "NITER", n, "actual number of iterations" );
    head.append<std::string>( "XREASON", reason, "reason for convergence" );
    head.append<realT>( "FPMTRANS", trans, "transmission of FPM" );
 
    fits::fitsFile<double> ff;
 
    std::string fname = "coron/" + cname + "_apod.fits";
 
    ff.write( fname, m_pupilApodizer, head );
 
    fname = "coron/" + cname + "_fpm.fits";
    improc::eigenImage<fpmaskFloatT> fpm( m_focalMask.rows(), m_focalMask.cols() );
    for( int r = 0; r < m_focalMask.rows(); ++r )
    {
        for( int c = 0; c < m_focalMask.cols(); ++c )
        {
            fpm( r, c ) = m_focalMask( r, c ).real();
        }
    }
 
    ff.write( fname, fpm, head );
 
    fname = "coron/" + cname + "_lyot.fits";
    ff.write( fname, m_lyotStop, head );
}
 
} // namespace wfp
} // namespace mx
 
#endif //__lyotCoronagraph_hpp__