aoAtmosphere.hpp

Go to the documentation of this file.

/** \file aoAtmosphere.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Provides a class to specify atmosphere parameters.
 * \ingroup mxAO_files
 *
 */
 
#ifndef aoAtmosphere_hpp
#define aoAtmosphere_hpp
 
#include <iostream>
#include <numeric>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <algorithm>
 
#include "../../mxlib.hpp"
#include "../../mxError.hpp"
 
#include "aoConstants.hpp"
 
#include "../../math/constants.hpp"
 
#include "../../app/appConfigurator.hpp"
 
namespace mx
{
namespace AO
{
namespace analysis
{
 
/// A class to specify atmosphere parameters and perform related calculations.
/**
 * \todo layer outer scales need work
 * \todo layer PSD params isn't finished.  Need to push that to PSD itself, and manage things like wavelength
 * dependence.
 *
 * \tparam realT is the real floating type in which all calculations are performed.
 *
 * \ingroup mxAOAnalytic
 */
template <typename _realT>
class aoAtmosphere
{
  public:
    typedef _realT realT; ///< The real floating type in which all calculations are performed.
 
    /// Constructor
    aoAtmosphere();
 
  protected:
    realT m_r_0; ///< Fried's parameter, in m
 
    realT m_lam_0{ 0.5e-6 }; ///< Wavelength of Fried's parameter, in m
 
    std::vector<realT> m_layer_Cn2; ///< Vector of layer strengths.
 
    std::vector<realT> m_L_0; ///< The outer scale, in m
 
    std::vector<realT> m_l_0; ///< The inner scale of each layer, in m
 
    bool m_nonKolmogorov{ false }; ///< Flag indicating if non-Kolmogorov PSD parameters are used.
 
    std::vector<realT> m_beta{ 1 }; ///< The PSD normalization when in non-Kolmogorov mode.
 
    std::vector<realT> m_alpha{ 0 }; ///< The PSD exponent when in non-Kolmogorov mode.
 
    std::vector<realT> m_beta_0{ 0 }; ///< The PSD constant when in non-Kolmogorov mode.
 
    std::vector<realT> m_layer_z; ///< Vector of layer heights, in m, above the observatory.
 
    realT m_h_obs{ 0 }; ///< Height of the observatory above sea level, in m.
 
    realT m_H{ 8000 }; ///< The atmospheric scale height, in m.
 
    std::vector<realT> m_layer_v_wind; ///< Vector of layer wind speeds, in m/s.
 
    std::vector<realT> m_layer_dir; ///< Vector of layer wind directions, in radians.
 
    bool m_v_wind_updated{ false }; ///< whether or not m_v_wind has been updated after changes
 
    realT m_v_wind; ///< \f$ C_n^2 \f$ averaged windspeed
 
    realT m_dir_wind; ///< \f$ C_n^2 \f$ averaged direction
 
    bool m_z_mean_updated{ false }; ///< whether or not m_z_mean has been updated after changes
 
    realT m_z_mean; ///< \f$ C_n^2 \f$ averaged layer height
 
    /// Checks if layer vectors have consistent length.
    /**
     * \returns 0 if all layer vectors are the same length
     * \returns -1 if not, and prints an error.
     */
    int checkLayers();
 
  public:
    /** \name PSD Parameters
     * @{
     */
 
    /// Get the value of Fried's parameter r_0 at the reference wavelength lam_0.
    /**
     * \returns the curret value of m_r_0, in m.
     */
    realT r_0();
 
    /// Get the value of Fried's parameter r_0 at the specified wavelength.
    /**
     * \note This does not set the value of r_0.
     *
     * \returns the current value of m_r_0 * pow(lam, 6/5), in m.
     */
    realT r_0( const realT &lam /**< [in] the wavelength, in m, at which to calculate r_0*/ );
 
    /// Set the value of Fried's parameter and the reference wavelength.
    /** If the provided reference wavelength is <=0, then 0.5 microns is used.
     *
     */
    void r_0( const realT &r0, ///< [in] is the new value of r_0, m
              const realT &l0  ///< [in] is the new value of lam_0 (in m), if 0 then 0.5e-6 is the default.
    );
 
    /// Get the current value of the reference wavelength.
    /** This is the wavelength at which r_0 is specified.
     *
     * \returns the current value of m_lam_0, in m.
     */
    realT lam_0();
 
    /// Get the strength of a single layer.
    /**
     * \returns the value of m_layer_Cn2[n].
     */
    realT layer_Cn2( const int n /**< [in] specifies the layer. */ );
 
    /// Get the vector of layer strengths.
    /**
     * \returns a copy of the vector of layer strengths: m_layer_Cn2.
     */
    std::vector<realT> layer_Cn2();
 
    /// Set the vector layer strengths, possibly calculating r_0.
    /**
     * If a reference wavelength is specified (l0 > 0), then r_0 is set from the layer strengths
     * according to
     *
     * Regardless of what units the strengths are specified in, they are stored normalized so that
     * \f$ \sum_n C_n^2 = 1 \f$.
     *
     */
    void layer_Cn2( const std::vector<realT> &cn2, ///<  [in] is a vector containing the layer strengths
                    const realT l0 = 0 ///< [in] [optional] if l0 > 0, then r_0 is set from the layer strengths.
    );
 
    /// Get the value of the outer scale for a single layer.
    /**
     * \returns the current value of m_L_0[n], in m.
     */
    realT L_0( const size_t &n );
 
    /// Set the vector of layer outer scales.
    /**
     */
    void L_0( const std::vector<realT> &L0 /**< [in] is the new vector of outer scales, in m  */ );
 
    /// Get the vector of outer scales.
    /**
     * \returns a copy of the vector of layer outer scales, in m
     */
    std::vector<realT> L_0();
 
    /// Get the value of the inner scale for a single layer.
    /**
     * \returns the current value of m_l_0[n], in m.
     */
    realT l_0( const size_t &n );
 
    /// Set the vector of layer inner scales.
    /**
     */
    void l_0( const std::vector<realT> &l0 /**< [in] is the new vector of inner scales, in m*/ );
 
    /// Get the vector of inner scales.
    /**
     * \returns a copy of the vector of layer inner scales, in m
     */
    std::vector<realT> l_0();
 
    /// Set the value of m_nonKolmogorov
    /** This flag indicates if non-Kolmogorov turbulence is being modeled.
     */
    void nonKolmogorov( const bool &nk /**< [in] the value of m_nonKolmogorov*/ );
 
    /// Return the value of m_nonKolmogorov
    /** This flag indicates if non-Kolmogorov turbulence is being modeled.
     *
     * \returns the current value of m_nonKolmogorov
     */
    bool nonKolmogorov();
 
    /// Return the PSD index for a single layer.
    /** Satifies the requirements of psdParamsT.
      *
      * If m_nonKolmogorov is false, this returns
      * \[
        \alpha = \frac{11}{3}
        \]
      * Otherwise it returns the current value of m_alpha[n].
      *
      * \returns the PSD index.
      */
    realT alpha( const size_t &n );
 
    /// Set the vector of layer PSD indices.
    /**
     */
    void alpha( const std::vector<realT> &alph /**< [in] is the new vector of PSD indices*/ );
 
    /// Get the vector of PSD indices.
    /**
     * \returns a copy of the vector of PSD indices
     */
    std::vector<realT> alpha();
 
    /// Return the PSD normalization for a single layer.
    /** Satifies the requirements of psdParamsT.
      *
      * If m_nonKolmogorov is false, this returns
      * \[
        \beta = \frac{0.0218}{r_0^{5/3}}
        \]
      * Otherwise it returns the current value of m_beta[n].
      *
      * \returns the PSD normalization.
      */
    realT beta( const size_t &n );
 
    /// Set the vector of layer PSD normalizations.
    /**
     */
    void beta( const std::vector<realT> &bet /**< [in] is the new vector of PSD normalizations*/ );
 
    /// Get the vector of PSD normalizations.
    /**
     * \returns a copy of the vector of PSD normalizations
     */
    std::vector<realT> beta();
 
    /// Return the PSD constant for a single layer.
    /** Satifies the requirements of psdParamsT.
      *
      * If m_nonKolmogorov is false, this returns
      * \[
        \beta_0 = 0
        \]
      * Otherwise it returns the current value of m_beta_0[n].
      *
      * \returns the PSD constant.
      */
    realT beta_0( const size_t &n );
 
    /// Set the vector of layer PSD constants.
    /**
     */
    void beta_0( const std::vector<realT> &bet /**< [in] is the new vector of PSD constants*/ );
 
    /// Get the vector of PSD constants.
    /**
     * \returns a copy of the vector of PSD constants
     */
    std::vector<realT> beta_0();
 
    /// Get the number of layers
    /**
     * \returns the size of the m_layer_Cn2 vector if m_nonKolmogrov==false
     * \returns the size of the m_alpha vector if m_nonKolmogrov==true
     */
    size_t n_layers();
 
    size_t currentLayer();
 
    void currentLayer( size_t cl );
 
    /// @}
 
    /// Get the height of a single layer.
    /**
     * \returns the height of layer n, in m.
     */
    realT layer_z( const size_t n /**< [in] specifies the layer. */ );
 
    /// Get the vector layer heights.
    /**
     * \returns a copy of the vector of layer heights:m_layer_z, in m.
     */
    std::vector<realT> layer_z();
 
    /// Set the vector of layer heights.
    /**
     */
    void layer_z( const std::vector<realT> &layz /**< [in] new vector of layer heights, in m. */ );
 
    /// Get the height of the observatory.
    /**
     * \return the current value of m_h_obs, in m.
     */
    realT h_obs();
 
    /// Set the height of the observatory.
    /**
     */
    void h_obs( realT nh /**< [in] the new height of the observatory [m] */ );
 
    /// Get the atmospheric scale height.
    /**
     * \returns the current value of m_H, in m.
     */
    realT H();
 
    /// Set the atmospheric scale height.
    /**
     */
    void H( realT nH /**< [in] the new value of m_H [m] */ );
 
    /// Get the wind speed of a single layer.
    /**
     * \returns the value of m_layer_v_wind[n].
     */
    realT layer_v_wind( const int n /**< [in] specifies the layer. */ );
 
    /// Get the vector of layer windspeeds
    /**
     * \returns a copy of the vector of layer windspeeds: m_layer_v_wind.
     */
    std::vector<realT> layer_v_wind();
 
    /// Set the vector of layer windspeeds.
    /**
     * \param
     */
    void layer_v_wind( const std::vector<realT> &spd /**< [in] the new vector, which is copied to m_layer_v_wind. */ );
 
    /// Get the wind direction of a single layer.
    /**
     * \returns the value of m_layer_dir[n], which is the wind direction in that layer in radians.
     */
    realT layer_dir( const int n /**< [in] specifies the layer. */ );
 
    /// Get the vector of layer wind directions
    /**
     * \returns a copy of the vector of layer wind directions: m_layer_dir.
     */
    std::vector<realT> layer_dir();
 
    /// Set the vector of layer wind directions.
    /**
     * \param
     */
    void layer_dir( const std::vector<realT>
                        &d /**< [in] the new vector of wind directions in radians, which is copied to m_layer_dir. */ );
 
    /// Get the 5/3 moment weighted mean wind speed
    /** Returns the weighted mean wind speed according to the 5/3's turbulence moment.  This is defined as
      *
      \f[
       \bar{v} = \left[\sum_i C_N^2(z_i) v_i^{5/3} \right]^{3/5}
      \f]
      * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
      *
      * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
      * returns the value of m_v_wind.
      *
      * \returns the current value of m_v_wind.
      */
    realT v_wind();
 
    /// Get the mean wind speed
    /** Returns the weighted mean wind speed according to the 5/3's turbulence moment.  This is defined as
      *
      \f[
       \bar{v} = \sum_i C_N^2(z_i) v_i
      \f]
      * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
      *
      * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
      * returns the value of m_v_wind.
      *
      * \returns the current value of m_v_wind.
      */
    realT v_wind_mean();
 
    /// Get the mean-squared wind speed
    realT v_wind_mean2();
 
    realT v_max()
    {
        auto max = std::max_element( std::begin( m_layer_v_wind ), std::end( m_layer_v_wind ) );
 
        return *max;
    }
 
    /// Get the weighted mean wind direction
    /** Returns the weighted mean wind direction according to the 5/3's turbulence moment.  This is defined as
      *
      \f[
       \bar{\theta} = \left[\sum_i C_N^2(z_i) \theta_i^{5/3} \right]^{3/5}
      \f]
      * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
      *
      * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
      * returns the value of m_dir_wind.
      *
      * \returns the current value of m_dir_wind.
      */
    realT dir_wind();
 
  protected:
    /// Recalculate m_v_wind
    /** Called by v_wind() whenever m_v_wind_updated is true.
     *
     * \todo handle dir_wind averaging across 0/2pi
     */
    void update_v_wind();
 
  public:
    /// Set the weighted mean m_v_wind and renormalize the layer wind speeds
    /** Calling this function changes the values of m_layer_v_wind so that the
     * Layer averaged 5/3 \f$C_n^2\f$ moment of wind speed is the new value specified by vw.
     *
     */
    void v_wind( const realT &vw /**< [in] the new value of m_v_wind. */ );
 
    /// Get the weighted mean layer height
    /** Returns the weighted layer height according to the 5/3's turbulence moment.  This is defined as
      \f[
       \bar{z} = \left[\sum_i C_N^2(z_i) z_i^{5/3} \right]^{3/5}
      \f]
      * See Hardy (1998) Section 3.3.6 and 3.7.2. \cite hardy_1998.
      *
      * This is only re-calculated if either m_layer_Cn2 orm_layer_z is changed, otherwise this just
      * returns the value of m_z_mean.
      *
      * \returns the current value of m_z_mean.
      */
    realT z_mean();
 
  protected:
    /// Recalculate m_z_mean
    /** Called by z_mean() whenever m_z_mean_updated is true.
     */
    void update_z_mean();
 
  public:
    /// Set the weighted mean m_z_mean and renormalize the layer heights
    /** Calling this function changes the values ofm_layer_z so that the
     * Layer averaged 5/3 \f$C_n^2\f$ moment of the height is the new value specified.
     *
     */
    void z_mean( const realT &zm /**< [in] is the new value of m_v_wind. */ );
 
    /// The fraction of the turbulence PSD in phase after Fresnel propagation.
    /** See Equation (14) of Guyon (2005) \cite guyon_2005.
     *
     * \returns the value of the X function.
     */
    realT X( realT k,       ///< [in] the spatial frequency, in inverse meters.
             realT lam_sci, ///< [in] is the science observation wavelength.
             realT secZ     ///< [in] is the secant of the zenith distance.
    );
 
    /// The differential fraction of the turbulence PSD in phase after Fresnel propagation.
    /** See Equation (25) of Guyon (2005) \cite guyon_2005.
     *
     * \returns the value of the dX function.
     */
    realT dX( realT k,       ///< [in] the spatial frequency, in inverse meters.
              realT lam_sci, ///<  [in] is the science observation wavelength.
              realT lam_wfs  ///< [in] is the wavefront sensor wavelength.
    );
 
    /// The fraction of the turbulence PSD in amplitude after Fresnel propagation.
    /** See Equation (15) of Guyon (2005) \cite guyon_2005.
     *
     * \returns the value of the Y function.
     */
    realT Y( realT k,       ///< [in] the spatial frequency, in inverse meters.
             realT lam_sci, ///< [in] is the science observation wavelength.
             realT secZ     ///< [in] is the secant of the zenith distance.
    );
 
    /// The differential fraction of the turbulence PSD in amplitude after Fresnel propagation.
    /** See Equation (27) of Guyon (2005) \cite guyon_2005.
     *
     * \returns the value of the dY function.
     */
    realT dY( realT k,       ///< [in] the spatial frequency, in inverse meters.
              realT lam_sci, ///< [in] is the science observation wavelength.
              realT lam_wfs  ///< [in] is the wavefront sensor wavelength.
    );
 
    realT n_air( realT lam /**< [in]The wavelength*/ );
 
    realT X_Z( realT k,          ///< [in] the spatial frequency, in inverse meters
               realT lambda_sci, ///< [in] is the science observation wavelength.
               realT lambda_wfs, ///< [in] is the wavefront sensor wavelength.
               realT secZ        ///< [in] is the secant of the zenith distance.
    );
 
    /// Calculate the full-width at half-maximum of a seeing limited image for this atmosphere for a small telescope
    /// (ignoring L_0)
    /** Calculate the FWHM of a seeing limited image with the current parameters according to Floyd et al. (2010) \cite
      floyd_2010 \f[ \epsilon_0 = 0.98\frac{\lambda_{sci}}{r_0(\lambda_sci)}. \f]
      *
      * \returns the value of the FWHM for the current atmosphere parameters.
      */
    realT fwhm0( realT lam_sci /**< [in] the wavelength of the science observation. */ );
 
    /// Calculate the full-width at half-maximum of a seeing limited image for this atmosphere for a large telescope
    /// (including L_0)
    /** Calculate the FWHM of a seeing limited image with the current parameters according to Floyd et al. (2010) \cite
      floyd_2010 \f[ \epsilon_0 = 0.98\frac{\lambda_{sci}}{r_0(\lambda_sci)}. \f]
      * If there is an outer scale (_L_0 > 0), then a correction is applied according to Tokovinin (2002) \cite
      tokovinin_2002 \f[ \left( \frac{\epsilon_{vK}}{\epsilon_0}\right)^2 = 1 - 2.183\left(
      \frac{r_0(\lambda_{sci}}{L_0}\right)^{0.356} \f]
      *
      *
      * \returns the value of the FWHM (\f$ \epsilon_{0/vK} \f$) for the current atmosphere parameters.
      */
    realT fwhm( realT lam_sci /**< [in] the wavelength of the science observation. */ );
 
    /// Get the greenwood frequency at the reference wavelength
    /**
     * \todo derive full value of the constant
     */
    realT f_g();
 
    /// Get the greenwood frequency at a specified wavelength
    /**
     *
     */
    realT f_g( realT lam_sci /**< [in] the wavelength of the science observation. */ );
 
    /// Get tau_0 at the reference wavelength
    /**
     * \todo derive full value of the constant
     */
    realT tau_0();
 
    /// Get tau_0 at a specified wavelength.
    /**
     *
     */
    realT tau_0( realT lam_sci /**< [in] the wavelength of the science observation. */ );
 
    /// Scale v_wind so that tau_0 has the specified value at the specified wavelength.
    /** Does not modify r_0.
     */
    void tau_0( realT tau_0,  ///< [in] the desired tau_0
                realT lam_sci ///< [in] the wavelength of the science observation.
    );
 
    /// Load the default atmosphere model from Guyon (2005).
    /** Sets the parameters from Table 4 of Guyon (2005) \cite guyon_2005.
     */
    void loadGuyon2005();
 
    /// Load parameters corresponding to the median atmosphere of the GMT site survey at LCO.
    /**
     */
    void loadLCO();
 
    /// Set a single layer model.
    /** Sets all layer vectors to size=1 and populates their fields based on these arguments.
     *
     */
    void setSingleLayer( realT r0,   ///< [in] is the new value of r_0
                         realT lam0, ///< [in] is the new value of lam_0, if 0 then 0.5 microns is the default.
                         realT L0,   ///< [in] the new outer scale
                         realT l0,   ///< [in] the new inner scale
                         realT lz,   ///< [in] the layer height
                         realT vw,   ///< [in] the layer wind-speed
                         realT dir   ///< [in] the layer wind direction.
    );
 
    /// Output current parameters to a stream
    /** Prints a formatted list of all current parameters.
     *
     * \tparam iosT is a std::ostream-like type.
     *
     * \todo update for new vector components (L_0, etc.)
     */
    template <typename iosT>
    iosT &dumpAtmosphere( iosT &ios /**< [in] a std::ostream-like stream. */ );
 
    /// \name mx::application support
    /** @{
     */
 
    /// Setup the configurator to configure this class
    /**
     * Tests:
     *     - Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]"
     */
    void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
 
    /// Load the configuration of this class from a configurator
    /**
     * Tests:
     *     - Loading aoAtmosphere config settings \ref tests_ao_analysis_aoAtmosphere_config "[test doc]"
     */
    void loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
 
    /// @}
};
 
template <typename realT>
aoAtmosphere<realT>::aoAtmosphere()
{
}
 
template <typename realT>
int aoAtmosphere<realT>::checkLayers()
{
    size_t n = m_L_0.size();
 
    if( m_l_0.size() != n )
    {
        mxError( "aoAtmosphere", MXE_SIZEERR, "mismatched layer numbers (inner scale vs. outer scale)" );
        return -1;
    }
 
    if( m_layer_z.size() != n )
    {
        mxError( "aoAtmosphere", MXE_SIZEERR, "mismatched layer numbers (layer_z  vs. outer scale)" );
        return -1;
    }
 
    if( m_layer_Cn2.size() != n )
    {
        mxError( "aoAtmosphere", MXE_SIZEERR, "mismatched layer numbers (layer_Cn2  vs. outer scale)" );
        return -1;
    }
 
    if( m_layer_dir.size() != n )
    {
        mxError( "aoAtmosphere", MXE_SIZEERR, "mismatched layer numbers (layer_dir  vs. outer scale)" );
        return -1;
    }
 
    if( m_layer_v_wind.size() != n )
    {
        mxError( "aoAtmosphere", MXE_SIZEERR, "mismatched layer numbers (layer_v_wind vs. outer scale)" );
        return -1;
    }
 
    return 0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::r_0()
{
    return m_r_0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::r_0( const realT &lam )
{
    return m_r_0 * pow( lam / m_lam_0, math::six_fifths<realT>() );
}
 
template <typename realT>
void aoAtmosphere<realT>::r_0( const realT &r0, const realT &l0 )
{
    m_r_0 = r0;
 
    if( l0 > 0 )
    {
        m_lam_0 = l0;
    }
    else
    {
        m_lam_0 = 0.5e-6;
    }
}
 
template <typename realT>
realT aoAtmosphere<realT>::lam_0()
{
    return m_lam_0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::layer_Cn2( const int n )
{
    return m_layer_Cn2[n];
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::layer_Cn2()
{
    return m_layer_Cn2;
}
 
template <typename realT>
void aoAtmosphere<realT>::layer_Cn2( const std::vector<realT> &cn2, const realT l0 )
{
    m_layer_Cn2 = cn2;
 
    realT layer_norm = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        layer_norm += cn2[i];
    }
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
        m_layer_Cn2[i] = m_layer_Cn2[i] / layer_norm;
 
    if( l0 > 0 )
    {
        m_r_0 = 1.0 / pow( layer_norm * 5.520e13, math::three_fifths<realT>() );
        m_lam_0 = l0;
    }
 
    m_v_wind_updated = false;
    m_z_mean_updated = false;
}
 
template <typename realT>
realT aoAtmosphere<realT>::L_0( const size_t &n )
{
    return m_L_0[n];
}
 
template <typename realT>
void aoAtmosphere<realT>::L_0( const std::vector<realT> &L_0 )
{
    m_L_0 = L_0;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::L_0()
{
    return m_L_0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::l_0( const size_t &n )
{
    return m_l_0[n];
}
 
template <typename realT>
void aoAtmosphere<realT>::l_0( const std::vector<realT> &l_0 )
{
    m_l_0 = l_0;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::l_0()
{
    return m_l_0;
}
 
template <typename realT>
void aoAtmosphere<realT>::nonKolmogorov( const bool &nk )
{
    m_nonKolmogorov = nk;
}
 
template <typename realT>
bool aoAtmosphere<realT>::nonKolmogorov()
{
    return m_nonKolmogorov;
}
 
template <typename realT>
realT aoAtmosphere<realT>::alpha( const size_t &n )
{
    if( !m_nonKolmogorov )
    {
        return math::eleven_thirds<realT>();
    }
    else
    {
        return m_alpha[n];
    }
}
 
template <typename realT>
void aoAtmosphere<realT>::alpha( const std::vector<realT> &alph )
{
    m_nonKolmogorov = true;
    m_alpha = alph;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::alpha()
{
    return m_alpha;
}
 
template <typename realT>
realT aoAtmosphere<realT>::beta( const size_t &n )
{
    if( !m_nonKolmogorov )
    {
        return constants::a_PSD<realT>() * pow( m_r_0, -math::five_thirds<realT>() );
        ;
    }
    else
    {
        return m_beta[n];
    }
}
 
template <typename realT>
void aoAtmosphere<realT>::beta( const std::vector<realT> &bet )
{
    m_nonKolmogorov = true;
    m_beta = bet;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::beta()
{
    return m_beta;
}
 
template <typename realT>
realT aoAtmosphere<realT>::beta_0( const size_t &n )
{
    if( !m_nonKolmogorov )
    {
        return 0;
    }
    else
    {
        return m_beta_0[n];
    }
}
 
template <typename realT>
void aoAtmosphere<realT>::beta_0( const std::vector<realT> &bet )
{
    m_nonKolmogorov = true;
    m_beta_0 = bet;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::beta_0()
{
    return m_beta_0;
}
 
template <typename realT>
size_t aoAtmosphere<realT>::n_layers()
{
    if( m_nonKolmogorov )
        return m_alpha.size();
    else
        return m_layer_Cn2.size();
}
 
template <typename realT>
realT aoAtmosphere<realT>::layer_z( const size_t n )
{
    return m_layer_z[n];
}
 
template <typename realT>
void aoAtmosphere<realT>::layer_z( const std::vector<realT> &z )
{
    m_layer_z = z;
    m_z_mean_updated = false;
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::layer_z()
{
    return m_layer_z;
}
 
template <typename realT>
realT aoAtmosphere<realT>::h_obs()
{
    return m_h_obs;
}
 
template <typename realT>
void aoAtmosphere<realT>::h_obs( realT nh )
{
    m_h_obs = nh;
}
 
template <typename realT>
realT aoAtmosphere<realT>::H()
{
    return m_H;
}
 
template <typename realT>
void aoAtmosphere<realT>::H( realT nH )
{
    m_H = nH;
}
 
template <typename realT>
realT aoAtmosphere<realT>::layer_v_wind( const int n )
{
    return m_layer_v_wind[n];
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::layer_v_wind()
{
    return m_layer_v_wind;
}
 
template <typename realT>
void aoAtmosphere<realT>::layer_v_wind( const std::vector<realT> &v )
{
    m_layer_v_wind = v;
    m_v_wind_updated = false;
}
 
template <typename realT>
realT aoAtmosphere<realT>::layer_dir( const int n )
{
    return m_layer_dir[n];
}
 
template <typename realT>
std::vector<realT> aoAtmosphere<realT>::layer_dir()
{
    return m_layer_dir;
}
 
template <typename realT>
void aoAtmosphere<realT>::layer_dir( const std::vector<realT> &dir )
{
    m_layer_dir = dir;
    m_v_wind_updated = false;
}
 
template <typename realT>
realT aoAtmosphere<realT>::v_wind()
{
    if( m_v_wind_updated == false )
        update_v_wind();
    return m_v_wind;
}
 
template <typename realT>
realT aoAtmosphere<realT>::dir_wind()
{
    if( m_v_wind_updated == false )
        update_v_wind();
    return m_dir_wind;
}
 
template <typename realT>
void aoAtmosphere<realT>::update_v_wind()
{
 
    if( m_layer_v_wind.size() == 0 || m_layer_Cn2.size() == 0 )
    {
        m_v_wind = 0;
        m_dir_wind = 0;
        m_v_wind_updated = true;
        return;
    }
 
    m_v_wind = 0;
 
    realT s = 0;
    realT c = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        m_v_wind += m_layer_Cn2[i] * pow( m_layer_v_wind[i], math::five_thirds<realT>() );
        s += pow( m_layer_v_wind[i], math::five_thirds<realT>() ) *
             sin( m_layer_dir[i] ); // pow( sin(m_layer_dir[i]), math::five_thirds<realT>() );
        c += pow( m_layer_v_wind[i], math::five_thirds<realT>() ) *
             cos( m_layer_dir[i] ); // pow( cos(m_layer_dir[i]), math::five_thirds<realT>() );
    }
 
    m_v_wind = pow( m_v_wind, math::three_fifths<realT>() );
 
    // m_dir_wind = atan(pow(s, math::three_fifths<realT>()) / pow(c, math::three_fifths<realT>()));
    m_dir_wind = atan( s / c );
    if( m_dir_wind < 0 )
        m_dir_wind += math::pi<realT>();
 
    m_v_wind_updated = true;
}
 
template <typename realT>
void aoAtmosphere<realT>::v_wind( const realT &vw )
{
    if( m_v_wind_updated == false )
        update_v_wind();
 
    realT vw_old = m_v_wind;
 
    m_v_wind = vw;
 
    // Now update the layers if needed
    if( m_layer_v_wind.size() > 0 )
    {
        for( size_t i = 0; i < m_layer_v_wind.size(); ++i )
        {
            m_layer_v_wind[i] = m_layer_v_wind[i] * ( vw / vw_old );
        }
    }
}
 
template <typename realT>
realT aoAtmosphere<realT>::z_mean()
{
    if( m_z_mean_updated == false )
        update_z_mean();
    return m_z_mean;
}
 
template <typename realT>
void aoAtmosphere<realT>::update_z_mean()
{
    if( m_layer_z.size() == 0 || m_layer_Cn2.size() == 0 )
    {
        m_z_mean = 0;
        m_z_mean_updated = true;
        return;
    }
 
    m_z_mean = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        m_z_mean += m_layer_Cn2[i] * pow( m_layer_z[i], math::five_thirds<realT>() );
    }
 
    m_z_mean = pow( m_z_mean, math::three_fifths<realT>() );
 
    m_z_mean_updated = true;
}
 
template <typename realT>
void aoAtmosphere<realT>::z_mean( const realT &zm )
{
    if( m_z_mean_updated == false )
        update_z_mean();
 
    realT zh_old = m_z_mean;
 
    m_z_mean = zm;
 
    // Now update the layers if needed
    if( m_layer_z.size() > 0 )
    {
        for( size_t i = 0; i < m_layer_z.size(); ++i )
        {
            m_layer_z[i] = m_layer_z[i] * ( zm / zh_old );
        }
    }
}
 
template <typename realT>
realT aoAtmosphere<realT>::X( realT k, realT lam_sci, realT secZ )
{
    realT c = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        c += m_layer_Cn2[i] * pow( cos( math::pi<realT>() * k * k * lam_sci * m_layer_z[i] * secZ ), 2 );
    }
 
    return c;
}
 
template <typename realT>
realT aoAtmosphere<realT>::dX( realT f, realT lam_sci, realT lam_wfs )
{
    realT c = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        c += m_layer_Cn2[i] * pow( ( cos( math::pi<realT>() * f * f * lam_sci * m_layer_z[i] ) -
                                     cos( math::pi<realT>() * f * f * lam_wfs * m_layer_z[i] ) ),
                                   2 );
    }
 
    return c;
}
 
template <typename realT>
realT aoAtmosphere<realT>::Y( realT k, realT lam_sci, realT secZ )
{
    realT c = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        c += m_layer_Cn2[i] * pow( sin( math::pi<realT>() * k * k * lam_sci * m_layer_z[i] * secZ ), 2 );
    }
    return c;
}
 
template <typename realT>
realT aoAtmosphere<realT>::dY( realT f, realT lam_sci, realT lam_wfs )
{
    realT c = 0;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        c += m_layer_Cn2[i] * pow( ( sin( math::pi<realT>() * f * f * lam_sci * m_layer_z[i] ) -
                                     sin( math::pi<realT>() * f * f * lam_wfs * m_layer_z[i] ) ),
                                   2 );
    }
 
    return c;
}
 
template <typename realT>
realT aoAtmosphere<realT>::n_air( realT lambda )
{
    realT ll2 = static_cast<realT>( 1 ) / pow( lambda / 1e-6, 2 );
 
    return 1.0 + 8.34213e-5 + 0.0240603 / ( 130.0 - ll2 ) + 0.00015997 / ( 38.9 - ll2 );
}
 
template <typename realT>
realT aoAtmosphere<realT>::X_Z( realT k, realT lambda_i, realT lambda_wfs, realT secZ )
{
    realT c = 0;
    realT sinZ = sqrt( 1.0 - pow( 1.0 / secZ, 2 ) );
    realT tanZ = sinZ * secZ;
    realT x0 = ( n_air( lambda_wfs ) - n_air( lambda_i ) ) * m_H * tanZ * secZ;
    realT x;
 
    for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
    {
        x = x0 * ( 1 - exp( ( m_layer_z[i] + m_h_obs ) / m_H ) );
        c += m_layer_Cn2[i] * pow( cos( math::pi<realT>() * k * k * lambda_i * m_layer_z[i] * secZ ), 2 ) *
             pow( sin( math::pi<realT>() * x * k * cos( 0. * 3.14 / 180. ) ), 2 );
    }
 
    return 4 * c;
}
 
template <typename realT>
realT aoAtmosphere<realT>::fwhm0( realT lam_sci )
{
    realT r0lam = r_0( lam_sci );
 
    realT fwhm = 0.98 * ( lam_sci / r0lam );
 
    return fwhm;
}
 
template <typename realT>
realT aoAtmosphere<realT>::fwhm( realT lam_sci )
{
    realT r0lam = r_0( lam_sci );
 
    realT fwhm = 0.98 * ( lam_sci / r0lam );
 
    ///\todo this needs to handle layers with different L_0
    if( L_0( 0 ) > 0 )
        fwhm *= sqrt( 1 - 2.183 * pow( r0lam / L_0( 0 ), 0.356 ) );
 
    return fwhm;
}
 
template <typename realT>
realT aoAtmosphere<realT>::f_g()
{
    return 0.428 * v_wind() / m_r_0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::f_g( realT lam_sci )
{
    return 0.428 * pow( m_lam_0 / lam_sci, math::six_fifths<realT>() ) * v_wind() / m_r_0;
}
 
template <typename realT>
realT aoAtmosphere<realT>::tau_0()
{
    return 0.134 / f_g();
}
 
template <typename realT>
realT aoAtmosphere<realT>::tau_0( realT lam_sci )
{
    return 0.134 / f_g( lam_sci );
}
 
template <typename realT>
void aoAtmosphere<realT>::tau_0( realT tau_0, realT lam_sci )
{
    realT vw = ( 0.134 / tau_0 ) / 0.428 * m_r_0 * pow( lam_sci / m_lam_0, math::six_fifths<realT>() );
    v_wind( vw );
}
 
template <typename realT>
void aoAtmosphere<realT>::loadGuyon2005()
{
    layer_Cn2( { 0.2283, 0.0883, 0.0666, 0.1458, 0.3350, 0.1350 } );
    layer_z( { 500, 1000, 2000, 4000, 8000, 16000 } );
    layer_v_wind( { 10., 10., 10., 10., 10., 10. } );
    layer_dir( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } );
 
    r_0( 0.2, 0.5e-6 );
 
    h_obs( 4200 );
}
 
template <typename realT>
void aoAtmosphere<realT>::loadLCO()
{
    layer_Cn2( { 0.42, 0.029, 0.062, 0.16, 0.11, 0.10, 0.12 } );
    layer_z( { 250., 500., 1000., 2000., 4000., 8000., 16000. } );
    layer_v_wind( { 10.0, 10.0, 20.0, 20.0, 25.0, 30.0, 25.0 } );
    layer_dir( { 1.05, 1.05, 1.31, 1.31, 1.75, 1.92, 1.75 } );
 
    r_0( 0.16, 0.5e-6 );
 
    L_0( { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 } );
    l_0( { 0.0, 0, 0, 0, 0, 0, 0 } );
 
    h_obs( 2400 );
}
 
template <typename realT>
void aoAtmosphere<realT>::setSingleLayer( realT r0, realT lam0, realT L0, realT l0, realT lz, realT vw, realT dir )
{
    r_0( r0, lam0 );
    L_0( std::vector<realT>( { L0 } ) );
    l_0( std::vector<realT>( { l0 } ) );
    layer_Cn2( std::vector<realT>( { 1 } ) );
    layer_z( std::vector<realT>( { lz } ) );
    layer_v_wind( std::vector<realT>( { vw } ) );
    layer_dir( std::vector<realT>( { dir } ) );
}
 
template <typename realT>
template <typename iosT>
iosT &aoAtmosphere<realT>::dumpAtmosphere( iosT &ios )
{
    ios << "# Atmosphere Parameters:\n";
    ios << "#    nonKolmogorov = " << std::boolalpha << nonKolmogorov() << '\n';
    ios << "#    n_layers = " << n_layers() << '\n';
 
    if( !m_nonKolmogorov )
    {
        ios << "#    r_0 = " << r_0() << '\n';
        ios << "#    lam_0 = " << lam_0() << '\n';
        ios << "#    tau_0 = " << tau_0( lam_0() ) << '\n';
        ios << "#    FWHM = " << fwhm( lam_0() ) << '\n';
        ios << "#    layer_Cn2 = ";
        for( size_t i = 0; i < n_layers() - 1; ++i )
            ios << layer_Cn2()[i] << ", ";
        ios << layer_Cn2()[n_layers() - 1] << '\n';
    }
    if( m_nonKolmogorov )
    {
        ios << "#    alpha = ";
        for( size_t i = 0; i < n_layers() - 1; ++i )
            ios << alpha()[i] << ", ";
        ios << alpha()[n_layers() - 1] << '\n';
        ios << "#    beta = ";
        for( size_t i = 0; i < n_layers() - 1; ++i )
            ios << beta()[i] << ", ";
        ios << beta()[n_layers() - 1] << '\n';
        ios << "#    beta_0 = ";
        for( size_t i = 0; i < n_layers() - 1; ++i )
            ios << beta_0()[i] << ", ";
        ios << beta_0()[n_layers() - 1] << '\n';
    }
 
    ios << "#    L_0 = ";
    for( size_t i = 0; i < n_layers() - 1; ++i )
        ios << L_0()[i] << ", ";
    ios << L_0()[n_layers() - 1] << '\n';
    ios << "#    l_0 = ";
    for( size_t i = 0; i < n_layers() - 1; ++i )
        ios << l_0()[i] << ", ";
    ios << l_0()[n_layers() - 1] << '\n';
 
    ios << "#    layer_v_wind = ";
    for( size_t i = 0; i < n_layers() - 1; ++i )
        ios << layer_v_wind()[i] << ", ";
    ios << layer_v_wind()[n_layers() - 1] << '\n';
    ios << "#    layer_dir = ";
    for( size_t i = 0; i < n_layers() - 1; ++i )
        ios << layer_dir()[i] << ", ";
    ios << layer_dir()[n_layers() - 1] << '\n';
    ios << "#    mean v_wind = " << v_wind() << '\n';
    ios << "#    mean dir_wind = " << dir_wind() << '\n';
 
    if( !m_nonKolmogorov )
    {
        ios << "#    layer_z = ";
        for( size_t i = 0; i < n_layers() - 1; ++i )
            ios << layer_z()[i] << ", ";
        ios << layer_z()[n_layers() - 1] << '\n';
        ios << "#    mean z = " << z_mean() << '\n';
        ios << "#    h_obs = " << h_obs() << '\n';
        ios << "#    H = " << H() << '\n';
    }
 
    return ios;
}
 
template <typename realT>
void aoAtmosphere<realT>::setupConfig( app::appConfigurator &config )
{
    using namespace mx::app;
 
    config.add( "atm.r_0", "", "atm.r_0", argType::Required, "atm", "r_0", false, "real", "Fried's parameter [m]" );
    config.add( "atm.lam_0",
                "",
                "atm.lam_0",
                argType::Required,
                "atm",
                "lam_0",
                false,
                "real",
                "The reference wavlength for r_0 [m]" );
    config.add(
        "atm.L_0", "", "atm.L_0", argType::Required, "atm", "L_0", false, "vector<real>", "Layer outer scales [m]" );
    config.add(
        "atm.l_0", "", "atm.l_0", argType::Required, "atm", "l_0", false, "vector<real>", "Layer inner scales [m]" );
    config.add( "atm.layer_z",
                "",
                "atm.layer_z",
                argType::Required,
                "atm",
                "layer_z",
                false,
                "vector<real>",
                "layer heights [m]" );
    config.add(
        "atm.h_obs", "", "atm.h_obs", argType::Required, "atm", "h_obs", false, "real", "height of observatory [m]" );
    config.add( "atm.H", "", "atm.H", argType::Required, "atm", "H", false, "real", "atmospheric scale heights [m]" );
    config.add( "atm.layer_Cn2",
                "",
                "atm.layer_Cn2",
                argType::Required,
                "atm",
                "layer_Cn2",
                false,
                "vector<real>",
                "Layer Cn^2.  Note that this does not set r_0." );
    config.add( "atm.layer_v_wind",
                "",
                "atm.layer_v_wind",
                argType::Required,
                "atm",
                "layer_v_wind",
                false,
                "vector<real>",
                "Layer wind speeds [m/s]" );
    config.add( "atm.layer_dir",
                "",
                "atm.layer_dir",
                argType::Required,
                "atm",
                "layer_dir",
                false,
                "vector<real>",
                "Layer wind directions [rad]" );
    config.add( "atm.v_wind",
                "",
                "atm.v_wind",
                argType::Required,
                "atm",
                "v_wind",
                false,
                "real",
                "Mean windspeed (5/3 momement), rescales layer speeds [m/s]" );
    config.add( "atm.tau_0",
                "",
                "atm.tau_0",
                argType::Required,
                "atm",
                "tau_0",
                false,
                "real",
                "AO time constant, sets v_wind and rescales layer speeds. [s]" );
    config.add( "atm.z_mean",
                "",
                "atm.z_mean",
                argType::Required,
                "atm",
                "z_mean",
                false,
                "real",
                "Mean layer height (5/3 momemnt), rescales layer heights [m/s]" );
    config.add( "atm.nonKolmogorov",
                "",
                "atm.nonKolmogorov",
                argType::Required,
                "atm",
                "nonKolmogorov",
                false,
                "bool",
                "Set to use a non-Kolmogorov PSD. See alpha and beta." );
    config.add( "atm.alpha",
                "",
                "atm.alpha",
                argType::Required,
                "atm",
                "alpha",
                false,
                "vector<real>",
                "Non-kolmogorov PSD exponent." );
    config.add( "atm.beta",
                "",
                "atm.beta",
                argType::Required,
                "atm",
                "beta",
                false,
                "vector<real>",
                "Non-kolmogorov PSD normalization." );
    config.add( "atm.beta_0",
                "",
                "atm.beta_0",
                argType::Required,
                "atm",
                "beta_0",
                false,
                "vector<real>",
                "Non-kolmogorov PSD constant." );
}
 
template <typename realT>
void aoAtmosphere<realT>::loadConfig( app::appConfigurator &config )
{
    // Here "has side effecs" means that the set function does more than simply copy the value.
 
    // The order of lam_0, Cn2, and r_0 is so that r_0 overrides the value set with Cn2 if lam_0 != 0.
    // lam_0 comes first because it calibrates r0 and Cn2
 
    // lam_0
    config( m_lam_0, "atm.lam_0" );
 
    // layer_Cn2
    std::vector<realT> lcn2 = m_layer_Cn2;
    config( lcn2, "atm.layer_Cn2" );
    if( config.isSet( "atm.layer_Cn2" ) )
        layer_Cn2( lcn2 );
 
    realT r0 = r_0();
    config( r0, "atm.r_0" );
    if( config.isSet( "atm.r_0" ) )
        r_0( r0, m_lam_0 );
 
    config( m_L_0, "atm.L_0" );
 
    config( m_l_0, "atm.l_0" );
 
    // Has side effects:
    std::vector<realT> layz = m_layer_z;
    config( layz, "atm.layer_z" ); // Do this no matter what to record source
    if( config.isSet( "atm.layer_z" ) )
        layer_z( layz ); // but only call this if changed
 
    config( m_h_obs, "atm.h_obs" );
    config( m_H, "atm.H" );
 
    // Has side effects:
    std::vector<realT> lvw = m_layer_v_wind;
    config( lvw, "atm.layer_v_wind" ); // Do this no matter what to record source
    if( config.isSet( "atm.layer_v_wind" ) )
        layer_v_wind( lvw ); // but only call this if changed
 
    // Has side effects:
    std::vector<realT> ld = m_layer_dir;
    config( ld, "atm.layer_dir" ); // Do this no matter what to record source
    if( config.isSet( "atm.layer_dir" ) )
        layer_dir( ld ); // but only call this if changed
 
    realT vw = m_v_wind;
    config( vw, "atm.v_wind" ); // Do this no matter what to record source
    if( config.isSet( "atm.v_wind" ) )
        v_wind( vw ); // but only call this if changed
 
    realT t0 = tau_0();
    config( t0, "atm.tau_0" ); // Do this no matter what to record source
    if( config.isSet( "atm.tau_0" ) )
    {
        std::cerr << "setting tau_0 " << t0 << "\n";
        tau_0( t0, m_lam_0 ); // but only call this if changed
    }
    realT zm = m_z_mean;
    config( zm, "atm.z_mean" ); // Do this no matter what to record source
    if( config.isSet( "atm.z_mean" ) )
        z_mean( zm ); // but only call this if changed
 
    config( m_nonKolmogorov, "atm.nonKolmogorov" );
 
    std::vector<realT> a = m_alpha;
    config( a, "atm.alpha" );
    if( config.isSet( "atm.alpha" ) )
        alpha( a ); // this sets m_nonKolmogorov
 
    std::vector<realT> b = m_beta;
    config( b, "atm.beta" );
    if( config.isSet( "atm.beta" ) )
        beta( b ); // this sets m_nonKolmogorov
 
    std::vector<realT> b0 = m_beta_0;
    config( b0, "atm.beta_0" );
    if( config.isSet( "atm.beta_0" ) )
        beta_0( b0 ); // this sets m_nonKolmogorov
}
 
extern template class aoAtmosphere<float>;
 
extern template class aoAtmosphere<double>;
 
extern template class aoAtmosphere<long double>;
 
#ifdef HASQUAD
extern template class aoAtmosphere<__float128>;
#endif
 
} // namespace analysis
} // namespace AO
} // namespace mx
 
#endif // aoAtmosphere_hpp