turbLayer.hpp

Go to the documentation of this file.

/** \file turbLayer.hpp
 * \brief Declaration and definition of a turbulence layer.
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup mxAO_sim_files
 *
 */
 
//***********************************************************************//
// Copyright 2023 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 mx_AO_sim_turbLayer_hpp
#define mx_AO_sim_turbLayer_hpp
 
#include <vector>
 
#include <Eigen/Dense>
 
#include "../../error/error.hpp"
#include "../../error/exception.hpp"
#include "../../improc/imageTransforms.hpp"
#include "../../math/constants.hpp"
#include "../../math/func/jinc.hpp"
#include "../../math/randomT.hpp"
#include "../../sigproc/psdFilter.hpp"
#include "../../sigproc/psdUtils.hpp"
 
namespace mx
{
namespace AO
{
namespace sim
{
 
template <typename _aoSystemT, class _verboseT>
struct turbAtmosphere;
 
/// Simulation of a single turbulent layer
/** \todo document this
 * \todo add facility for changing interpolator
 *
 * \ingroup mxAOSim
 */
template <typename _aoSystemT, class _verboseT = mx::verbose::d>
struct turbLayer
{
    typedef _aoSystemT aoSystemT;
    typedef _verboseT verboseT;
 
    typedef typename aoSystemT::realT realT;
    typedef Eigen::Array<realT, -1, -1> imageT;
 
    turbAtmosphere<aoSystemT, verboseT> *m_parent{ nullptr };
 
    uint32_t m_layerNo;
 
    uint32_t m_scrnSz;
 
    imageT m_freq;
    imageT m_psd;
 
    realT m_dx{ 0 };
    realT m_dy{ 0 };
 
    realT m_x0;
    realT m_y0;
 
    int m_last_wdx;
    int m_last_wdy;
 
    imageT m_phase;
    imageT m_shiftPhaseWP;
    imageT m_shiftPhase;
    imageT m_shiftPhaseWork;
 
    mx::math::uniDistT<realT> uniVar; ///< Uniform deviate, used in shiftRandom.
 
    void setLayer( turbAtmosphere<aoSystemT, verboseT> *parent, int layerNo, uint32_t scrnSz );
 
    uint32_t layerNo();
 
    uint32_t scrnSz();
 
    void alloc();
 
    /// Deallocate memory necessary for phase screen generation.
    void genDealloc();
 
    /// Shift to a timestep.
    /**
     * \param [in] dt is the new timestep.
     */
    void shift( realT dt );
 
    /// Seed the uniform deviation.  Call this if you intend to use shiftRandom.
    /** This only needs to be called once.
     */
    void initRandom();
 
    /// Shift by a random amount using the uniform distribution
    /** Call initRandom() once before calling this method.
     *
     * \param [in] nofract if true then the fractional part is ignored, and only a whole-pixel shift is executed.
     * Default=false
     */
    void shiftRandom( bool nofract = false );
};
 
// template<typename aoSystemT>
// turbLayer<aoSystemT, verboseT>::turbLayer()
// {
// }
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::setLayer( turbAtmosphere<aoSystemT, verboseT> *parent,
                                               int layerNo,
                                               uint32_t scrnSz )
{
    m_parent = parent;
    m_layerNo = layerNo;
 
    m_scrnSz = scrnSz;
 
    if( m_parent == nullptr )
    {
        throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent is nullptr)" );
    }
 
    if( m_parent->aosys() == nullptr )
    {
        throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent->aosys() is nullptr)" );
    }
 
    realT vwind = m_parent->aosys()->atm.layer_v_wind( m_layerNo );
    realT dwind = m_parent->aosys()->atm.layer_dir( m_layerNo );
    realT pupD = m_parent->aosys()->D();
    uint32_t wfSz = m_parent->wfSz();
 
    m_dx = vwind * cos( dwind ) / ( pupD / wfSz );
    m_dy = vwind * sin( dwind ) / ( pupD / wfSz );
 
    alloc();
}
 
template <typename aoSystemT, class verboseT>
uint32_t turbLayer<aoSystemT, verboseT>::layerNo()
{
    return m_layerNo;
}
 
template <typename aoSystemT, class verboseT>
uint32_t turbLayer<aoSystemT, verboseT>::scrnSz()
{
    return m_scrnSz;
}
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::alloc()
{
    if( m_parent == nullptr )
    {
        throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent is nullptr)" );
    }
 
    if( m_parent->aosys() == nullptr )
    {
        throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent->aosys() is nullptr)" );
    }
 
    m_phase.resize( m_scrnSz, m_scrnSz );
 
    uint32_t wfSz = m_parent->wfSz();
 
    uint32_t buffSz = m_parent->buffSz();
 
    m_shiftPhaseWP.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
 
    m_shiftPhase.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
 
    m_shiftPhaseWork.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
 
    initRandom();
 
    m_x0 = 0;
    m_y0 = 0;
    m_last_wdx = m_scrnSz + 1;
    m_last_wdy = m_scrnSz + 1;
 
    // Now set up for generation
    realT r0 = m_parent->aosys()->atm.r_0( m_parent->aosys()->lam_sci() );
    realT L0 = m_parent->aosys()->atm.L_0( m_layerNo );
    realT l0 = m_parent->aosys()->atm.l_0( m_layerNo );
 
    m_psd.resize( m_scrnSz, m_scrnSz );
 
    m_freq.resize( m_scrnSz, m_scrnSz );
    sigproc::frequencyGrid<imageT>( m_freq, m_parent->aosys()->D() / wfSz );
 
    realT beta = 0.0218 / pow( r0, 5. / 3. );
    realT sqrt_alpha = 0.5 * 11. / 3.;
 
    realT L02;
    if( L0 > 0 )
    {
        L02 = 1.0 / ( L0 * L0 );
    }
    else
    {
        L02 = 0;
    }
 
    // clang-format off
    #pragma omp parallel for // clang-format on
    for( size_t jj = 0; jj < m_scrnSz; ++jj )
    {
        for( size_t ii = 0; ii < m_scrnSz; ++ii )
        {
            realT p;
            if( m_freq( ii, jj ) == 0 && L02 == 0 )
            {
                p = 0;
            }
            else
            {
                p = beta / pow( pow( m_freq( ii, jj ), 2 ) + L02, sqrt_alpha );
                if( l0 > 0 )
                {
                    p *= exp( -1 * pow( m_freq( ii, jj ) * l0, 2 ) );
                }
            }
 
            realT Ppiston = 0;
            realT Ptiptilt = 0;
            if( m_parent->aosys()->psd.subPiston() )
            {
                Ppiston =
                    pow( 2 * math::func::jinc( math::pi<realT>() * m_freq( ii, jj ) * m_parent->aosys()->D() ), 2 );
            }
            if( m_parent->aosys()->psd.subTipTilt() )
            {
                Ptiptilt =
                    pow( 4 * math::func::jincN( 2, math::pi<realT>() * m_freq( ii, jj ) * m_parent->aosys()->D() ), 2 );
            }
 
            m_psd( ii, jj ) = sqrt( p * ( 1.0 - Ppiston - Ptiptilt ) );
        }
    }
 
    if( m_parent->shLevel() > 0 )
    {
        realT onoff = 1;
        if( m_parent->outerSubHarmonics() )
        {
            onoff = 0; // just make all of these 0
        }
 
        m_psd( 0, 0 ) = 0;
        m_psd( 0, 1 ) *= 0.5 * 0.5 * onoff;
        m_psd( 1, 0 ) *= 0.5 * 0.5 * onoff;
        m_psd( m_psd.rows() - 1, 0 ) *= 0.5 * 0.5 * onoff;
        m_psd( 0, m_psd.cols() - 1 ) *= 0.5 * 0.5 * onoff;
 
        m_psd( 1, 1 ) *= 0.75 * 0.75 * onoff;
        m_psd( m_psd.rows() - 1, m_psd.cols() - 1 ) *= 0.75 * 0.75 * onoff;
        m_psd( m_psd.rows() - 1, 1 ) *= 0.75 * 0.75 * onoff;
        m_psd( 1, m_psd.cols() - 1 ) *= 0.75 * 0.75 * onoff;
    }
}
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::genDealloc()
{
    m_freq.resize( 0, 0 );
    m_psd.resize( 0, 0 );
}
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::shift( realT dt )
{
    int wdx, wdy;
    realT ddx, ddy;
 
    ddx = m_x0 + m_dx * dt;
    ddy = m_y0 + m_dy * dt;
 
    wdx = static_cast<int>( std::floor( ddx ) );
    ddx -= wdx;
 
    wdy = static_cast<int>( std::floor( ddy ) );
    ddy -= wdy;
 
    wdx %= static_cast<int>( m_scrnSz );
    wdy %= static_cast<int>( m_scrnSz );
 
    if( wdx < 0 )
        wdx += m_scrnSz;
    if( wdy < 0 )
        wdy += m_scrnSz;
 
    if( dt == 0 )
    {
        m_shiftPhase = m_phase;
        // Prime both the whole-pixel and shifted buffers so the first
        // subsequent sub-pixel shift has valid input data.
        improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
        m_shiftPhase = m_shiftPhaseWP;
    }
    else
    {
        // Check for a new whole-pixel shift
        if( wdx != m_last_wdx || wdy != m_last_wdy || m_last_wdx == m_scrnSz + 1 || m_last_wdy == m_scrnSz + 1 )
        {
            // Need a whole pixel shift (this also extracts the m_wfSz + m_buffSz subarray)
            improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
        }
 
        // Do the sub-pixel shift
        improc::imageShift( m_shiftPhase, m_shiftPhaseWP, ddx, ddy, improc::cubicConvolTransform<realT>( -0.5 ) );
    }
 
    m_last_wdx = wdx;
    m_last_wdy = wdy;
}
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::initRandom()
{
    uniVar.seed();
}
 
template <typename aoSystemT, class verboseT>
void turbLayer<aoSystemT, verboseT>::shiftRandom( bool nofract )
{
    int wdx, wdy;
    realT ddx, ddy;
 
    ddx = uniVar * ( m_scrnSz );
    ddy = uniVar * ( m_scrnSz );
 
    wdx = (int)trunc( ddx );
    ddx -= wdx;
    wdy = (int)trunc( ddy );
    ddy -= wdy;
 
    if( nofract )
    {
        ddx = 0;
        ddy = 0;
    }
 
    improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
 
    if( ddx != 0 || ddy != 0 )
    {
        improc::imageShift( m_shiftPhase, m_shiftPhaseWP, ddx, ddy, improc::cubicConvolTransform<realT>( -0.5 ) );
    }
    else
    {
        m_shiftPhase = m_shiftPhaseWP;
    }
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif // mx_AO_sim_turbLayer_hpp