turbSequence.hpp

Go to the documentation of this file.

/** \file turbSequence.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Manage a turbulence sequence saved as FITS files.
 * \ingroup mxAO_sim_files
 *
 */
 
#ifndef turbSequence_hpp
#define turbSequence_hpp
 
#include <vector>
#include <string>
#include <iostream>
 
#include <Eigen/Dense>
 
#include "../../ioutils/fileUtils.hpp"
#include "../../ioutils/fits/fitsFile.hpp"
#include "../../error/mxErrorOld.hpp"
#include "../../improc/eigenCube.hpp"
 
#include "wavefront.hpp"
 
namespace mx
{
namespace AO
{
namespace sim
{
 
template <typename _realT>
struct turbSequence
{
    typedef _realT realT;
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
 
    std::vector<std::string> _phaseFnames;
    std::vector<std::string> _ampFnames;
 
    int _size{ 0 };
 
    int _nPerCube{ 0 };
 
    int _files{ 0 };
 
    int _frames{ 0 };
 
    int _currFileNo{ 0 };
    int _currFrameNo{ 0 };
 
    realT _wfPS{ 1 };
    realT _F0Photons{ 2e11 };
 
    realT _starMag{ 0 };
 
    realT _pixVal;
 
    imageT *_pupil;
 
    bool _loopClosed;
 
    improc::eigenCube<realT> _currPhase;
    improc::eigenCube<realT> _currAmp;
 
    bool _phaseOnly{ false };
 
    turbSequence()
    {
    }
 
    /// Return the size of the turbulence screens
    /**
     * \returns _size
     */
    int size()
    {
        return _size;
    }
 
    /// Return the number of frames per cube
    /**
     * \returns _nPerCube
     */
    int nPerCube()
    {
        return _nPerCube;
    }
 
    /// Return the number of files
    /**
     * \returns _files
     */
    int files()
    {
        return _files;
    }
 
    /// Return the number of frames in each file
    /**
     * \returns _frames
     */
    int frames()
    {
        return _frames;
    }
 
    realT wfPS()
    {
        return _wfPS;
    }
 
    void wfPS( realT ps )
    {
        _wfPS = ps;
        calcPixVal();
    }
 
    realT F0Photons()
    {
        return _F0Photons;
    }
 
    void F0Photons( realT f0 )
    {
        _F0Photons = f0;
        calcPixVal();
    }
 
    realT starMag()
    {
        return _starMag;
    }
 
    void starMag( realT sm )
    {
        _starMag = sm;
        calcPixVal();
    }
 
    realT pixVal()
    {
        return _pixVal;
    }
 
    void calcPixVal()
    {
        // This is the square root of photons/pixel
        _pixVal = sqrt( _F0Photons ) * pow( 10., -0.2 * _starMag ) * _wfPS; // pow(_wfPS,2);
    }
 
    /// Get the file names of the sequence assuming they are stored in a directory.
    /**
     * \param [in] dir is the directory to search for phase and amplitude files.
     * \param [in] max [optional] specifies the maximum number of files to include the file name list.
     *
     * \retval 0 on success.
     * \retval -1 on error.
     */
    int turbFnames( std::string dir, int max = 0 );
 
    void openPhaseFrame( int fn );
 
    void nextPhase( wavefront<realT> &wf );
 
    void nextWF( wavefront<realT> &wf );
};
 
template <typename realT>
int turbSequence<realT>::turbFnames( std::string dir, int max )
{
    ioutils::getFileNames( _phaseFnames, dir, "", ".pha", ".fits" );
    ioutils::getFileNames( _ampFnames, dir, "", ".amp", ".fits" );
 
    if( _phaseFnames.size() == 0 )
    {
        mxError( "turbSequence", MXE_FILENOTFOUND, "No turbulent phase files found." );
        return -1;
    }
 
    if( _ampFnames.size() == 0 )
    {
        std::cerr << "turbSequence: no turbulent amplitude files found.\n";
 
        _phaseOnly = true;
    }
 
    if( max > 0 )
    {
        const size_t maxFiles = static_cast<size_t>( max );
 
        if( _phaseFnames.size() > maxFiles )
            _phaseFnames.erase( _phaseFnames.begin() + maxFiles, _phaseFnames.end() );
 
        if( _phaseOnly == false && _ampFnames.size() > maxFiles )
            _ampFnames.erase( _ampFnames.begin() + maxFiles, _ampFnames.end() );
    }
 
    _files = _phaseFnames.size();
 
    openPhaseFrame( 0 );
 
    _size = _currPhase.rows();
    _nPerCube = _currPhase.planes();
 
    _frames = _nPerCube * _phaseFnames.size();
 
    return 0;
}
 
template <typename realT>
void turbSequence<realT>::openPhaseFrame( int fn )
{
    fits::fitsFile<float> ff;
 
    std::cout << _phaseFnames[fn] << "\n";
 
    ff.read( _currPhase, _phaseFnames[fn] );
    ff.close();
 
    if( _phaseOnly == false )
    {
        ff.read( _currAmp, _ampFnames[fn] );
        ff.close();
    }
 
    _currFileNo = fn;
    _currFrameNo = 0;
}
 
template <typename realT>
void turbSequence<realT>::nextPhase( wavefront<realT> &wf )
{
 
    int Npix = _pupil->sum();
 
    if( _currFrameNo > _nPerCube - 1 )
    {
        openPhaseFrame( _currFileNo + 1 );
    }
 
    wf.phase.resize( _pupil->rows(), _pupil->cols() );
    // mx::cutImage(wf.phase, _currPhase.image(_currFrameNo), _pupil->rows(), _pupil->cols());
 
    wf.phase = _currPhase.image( _currFrameNo )
                   .block( 0.5 * ( _currPhase.rows() - _pupil->rows() ),
                           0.5 * ( _currPhase.cols() - _pupil->cols() ),
                           _pupil->rows(),
                           _pupil->cols() );
 
    wf.phase = ( wf.phase - ( wf.phase * ( *_pupil ) ).sum() / Npix ) * ( *_pupil );
 
    // wf.amplitude.Constant(_pixVal);
 
    wf.amplitude = ( *_pupil ) * _pixVal;
 
    ++_currFrameNo;
}
 
template <typename realT>
void turbSequence<realT>::nextWF( wavefront<realT> &wf )
{
 
    int Npix = _pupil->sum();
 
    if( _currFrameNo > _nPerCube - 1 )
    {
        openPhaseFrame( _currFileNo + 1 );
    }
    wf.phase.resize( _pupil->rows(), _pupil->cols() );
    // mx::cutImage(wf.phase, _currPhase.image(_currFrameNo), _pupil->rows(), _pupil->cols());
    wf.phase = _currPhase.image( _currFrameNo )
                   .block( 0.5 * ( _currPhase.rows() - _pupil->rows() ),
                           0.5 * ( _currPhase.cols() - _pupil->cols() ),
                           _pupil->rows(),
                           _pupil->cols() );
 
    wf.phase = ( wf.phase - ( wf.phase * ( *_pupil ) ).sum() / Npix ) * ( *_pupil );
 
    if( _phaseOnly )
    {
        wf.amplitude = _pixVal * ( *_pupil );
    }
    else
    {
        wf.amplitude.resize( _pupil->rows(), _pupil->cols() );
        // mx::cutImage(wf.amplitude, _currAmp.image(_currFrameNo), _pupil->rows(), _pupil->cols());
        wf.amplitude = _currAmp.image( _currFrameNo )
                           .block( 0.5 * ( _currAmp.rows() - _pupil->rows() ),
                                   0.5 * ( _currAmp.cols() - _pupil->cols() ),
                                   _pupil->rows(),
                                   _pupil->cols() );
 
        wf.amplitude *= _pixVal * ( *_pupil );
    }
    _currFrameNo++;
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif // turbSequence_hpp