aoSystem.hpp

Go to the documentation of this file.

/** \file aoSystem.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Declares and defines an analytical AO system
 * \ingroup mxAO_files
 *
 */
 
//***********************************************************************//
// Copyright 2015-2022 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 aoSystem_hpp
#define aoSystem_hpp
 
#include "../../math/constants.hpp"
#include "../../math/roots.hpp"
#include "../../mxError.hpp"
#include "../../mxException.hpp"
 
#include "aoConstants.hpp"
 
#include "aoAtmosphere.hpp"
#include "aoPSDs.hpp"
#include "aoWFS.hpp"
 
#define FITTING_ERROR_NO 0
#define FITTING_ERROR_ZERO 1
#define FITTING_ERROR_X 2
#define FITTING_ERROR_Y 3
 
namespace mx
{
namespace AO
{
namespace analysis
{
 
/// Describes an analytic adaptive optics (AO) system.
/**
 * Templatized by the turbulence power spectral density (PSD).
 *
 * \tparam realT the floating point type used for all calculations
 * \tparam inputSpecT specifies the turbulence spatial PSD type
 * \tparam iosT is an output stream type with operator << defined (default is std::ostream)
 *
 * \ingroup mxAOAnalytic
 */
template <typename _realT, class _inputSpectT, typename iosT = std::ostream>
class aoSystem
{
 
  public:
    typedef _realT realT;
    typedef _inputSpectT inputSpectT;
    typedef aoAtmosphere<realT> aoAtmosphereT;
 
    aoAtmosphereT atm;
    inputSpectT psd;
 
  protected:
    realT m_D{ 0 }; ///< Telescope diameter [m]
 
    std::vector<realT> m_d_min{ 0 }; ///< Minimum AO system actuator pitch [m].  One per WFS mode.
 
    realT m_d_opt{ 0 };   ///< Current optimum AO system actuator pitch [m]
    bool m_optd{ false }; ///< Flag controlling whether actuator pitch is optimized (true) or just uses m_d_min (false).
                          ///< Default: true.
    realT m_optd_delta{ 1.0 }; ///< The fractional change from d_min used in optimization.  Set to 1 for integer
                               ///< binnings, > 1 for finer sampling.
 
    wfs<realT, iosT> *m_wfsBeta{ nullptr }; ///< The WFS beta_p class.
    bool m_ownWfsBeta{ false };             ///< Flag indicating if the WFS beta_p pointer is owned by this instance.
 
    realT m_opticalGain{
        1 }; ///< The optical gain of the WFS, 0-1. Treated as a sensitivity reduction in measurement noise.  Default 1.
 
  protected:
    realT m_lam_wfs{ 0 }; ///< WFS wavelength [m]
 
    // WFS detector parameters, which may be a function of binning
    std::vector<realT> m_npix_wfs{ 0 };  ///< Number of WFS pixels.  One per WFS mode.
    std::vector<realT> m_ron_wfs{ 0 };   ///< WFS readout noise [electrons/pix].  One per WFS mode.
    std::vector<realT> m_Fbg{ 0 };       ///< Background flux, [photons/sec/pixel].  One per WFS mode.
    std::vector<realT> m_minTauWFS{ 0 }; ///< Minimum WFS exposure time [sec].  One per WFS mode.
 
    bool m_bin_npix{ true }; ///< Flag controlling whether or not to bin WFS pixels according to the actuator spacing.
 
    int m_bin_opt{
        0 }; ///< The optimum binning factor.  If WFS modes are used, this is the mode index (0 to N-1).  If not, it is
             ///< 1 minus the pixel binning factor. It is always 1 minus the actuator binning factor.
 
    realT m_tauWFS{ 0 }; ///< Actual WFS exposure time [sec]
 
    realT m_deltaTau{ 0 }; ///< Loop latency [sec]
 
    bool m_optTau{ true }; ///< Flag controlling whether optimum integration time is calculated (true) enforcing
                           ///< m_minTauWFS, or if m_tauWFS is used (false). Default: true.
 
    realT m_lam_sci{ 0 }; ///< Science wavelength [m]
 
    realT m_zeta{ 0 };    ///<  Zenith angle [radians]
    realT m_secZeta{ 1 }; ///< Secant of the Zenith angle (calculated)
 
    int m_fit_mn_max{ 100 }; ///< Maximum spatial frequency index to use for fitting error calculation.
 
    bool m_circularLimit{ false }; ///< Flag to indicate that the spatial frequency limit is circular, not square.
 
    realT m_spatialFilter_ku{ std::numeric_limits<realT>::max() }; ///< The spatial filter cutoff in u, [m^-1]
    realT m_spatialFilter_kv{ std::numeric_limits<realT>::max() }; ///< The spatial filter cutoff in v, [m^-1]
 
    realT m_ncp_wfe{ 0 };   ///< Static WFE [m rms]
    realT m_ncp_alpha{ 2 }; ///< Power-law exponent for the NCP aberations.  Default is 2.
 
    realT m_F0{ 0 }; ///< 0 mag flux from star at WFS [photons/sec]
 
    realT m_starMag{ 0 }; ///< The magnitude of the star.
 
    bool m_specsChanged{ true }; ///< Flag to indicate that a specification has changed.
    bool m_dminChanged{ true };  ///< Flag to indicate that d_min has changed.
 
    Eigen::Array<int, -1, -1> m_controlledModes; ///< Map of which modes are under control.  Set by calcStrehl.
 
    realT m_wfeMeasurement{ 0 };   ///< Total WFE due to measurement a error [rad^2 at m_lam_sci]
    realT m_wfeTimeDelay{ 0 };     ///< Total WFE due to time delay [rad^2 at m_lam_sci]
    realT m_wfeFitting{ 0 };       ///< Total WFE due to fitting error [rad^2 at m_lam_sci]
    realT m_wfeChromScintOPD{ 0 }; ///< Total WFE due to the chromaticity of scintillation OPD [rad^2 at lam_sc]
    realT m_wfeChromIndex{ 0 };    ///< Total WFE due to the chromaticity of the index of refraction [rad^2 at lam_Sci]
    realT m_wfeAnisoOPD{ 0 };      ///< Total WFE due to dispersive anisoplanatism OPD.
 
    realT m_wfeNCP{ 0 }; ///< Total WFE due to NCP errors [rad^2 at m_lam_sci]
 
    realT m_wfeVar{ 0 }; ///< The WFE variance, in meters^2.  Never use this directly, instead use wfeVar().
    realT m_strehl{ 0 }; ///< Strehl ratio, a calculated quantity.  Never use this directdy, instead use strehl().
 
  public:
    /// Default c'tor
    /** Calls initialize().
     */
    aoSystem();
 
    /// Destructor
    ~aoSystem();
 
    /// Load the default parameters from Guyon, 2005 \cite guyon_2005.
    /**
     *
     */
    void loadGuyon2005();
 
    /// Load parameters corresponding to the MagAO-X system.
    void loadMagAOX();
 
    /// Load parameters corresponding to the G-MagAO-X system.
    void loadGMagAOX();
 
    /// Set the value of the primary mirror diameter.
    /**
     */
    void D( realT nD /**< [in] is the new value of m_D. */ );
 
    /// Get the value of the primary mirror diamter
    /**
     * \returns the current value of m_D.
     */
    realT D();
 
    /// Set the minimum subaperture sampling for each WFS mode.
    /** This sets the vector
     */
    void d_min( const std::vector<realT> &nd /**< [in] is the new values of m_d_min */ );
 
    /// Set the minimum subaperture sampling for a given WFS mode
    /** This sets the entry in the vector, if it has the proper length.
     */
    void d_min( int idx, ///< [in] the index of the WFS mode
                realT nd ///< [in] is the new value of m_d_min[idx]
    );
 
    /// Get the minimum subaperture sampling for a given WFS mode
    /**
     * \returns the current value of m_d_min[idx]
     */
    realT d_min( size_t idx /**< [in] the index of the WFS mode.*/ );
 
    /// Get the minimum subaperture sampling for all WFS modes.
    /**
     * \returns the current value of m_d_min.
     */
    std::vector<realT> d_min();
 
    /// Set whether or not the value of d is optimized or just set to m_d_min.
    /**
     */
    void optd( bool od /**< [in] is the new value of m_optd */ );
 
    /// Get the value of m_optd.
    /**
     * \returns the new value of m_optd_delta.
     */
    bool optd();
 
    /// Set the fractional change in actuator spacing for optimization.
    /** Sets the fraction of m_d_min by which the optimizer changes actautor spacing.
     */
    void optd_delta( realT odd /**< [in] is the new value of m_optd_delta */ );
 
    /// Get the value of the fractional change in actuator spacing for optimization..
    /**
     * \returns the value of m_optd_delta.
     */
    realT optd_delta();
 
    /// Set the WFS Beta pointer
    /** The WFS beta parameter determines the photon noise sensitivity.
     *
     * \tparam wfsT is a type derived from ao::analysis::wfs
     */
    template <typename wfsT>
    void wfsBeta( const wfsT &w /**< [in] an object derived from ao::analysis::wfs base class*/ );
 
    /// Set the WFS Beta pointer
    /** The WFS beta parameter determines the photon noise sensitivity.
     *
     * \tparam wfsT is a type derived from ao::analysis::wfs
     */
    template<typename wfsT>
   void wfsBeta( const wfsT * w /**< [in] pointer to an object derived from ao::analysis::wfs base class. If nullptr then a new object is allocated and managed.*/);
 
    /// Get the WFS Beta pointer
    /**
     * \returns a pointer to the current m_wfsBeta
     */
    wfs<realT, iosT> *wfsBeta();
 
    /// Get the value of beta_p for a spatial frequency
    /** beta_p is the photon noise sensitivity of the WFS.
     *
     * \returns beta_p as calculated by the WFS.
     */
    realT beta_p( realT m, ///< [in] the spatial frequency index
                  realT n  ///< [in] the spatial frequency index
    );
 
    /// Get the value of beta_r for a spatial frequency
    /** beta_r is the read noise sensitivity of the WFS.
     *
     * \returns beta_r as calculated by the WFS.
     */
    realT beta_r( realT m, ///< [in] the spatial frequency index
                  realT n  ///< [in] the spatial frequency index
    );
 
    /// Set the value of the Optical Gain.
    /**
     */
    void opticalGain( realT og /**< [in] the new value of m_opticalGain*/ );
 
    /// Get the value of the optical gain
    /** Optical gain is applied in the WFS noise calculation
     *
     * \returns the current value of m_opticalGain
     */
    realT opticalGain();
 
    /// Set the value of the WFS wavelength.
    /**
     */
    void lam_wfs( realT nlam /**< [in] is the new value of m_lam_wfs */ );
 
    /// Get the value of the WFS wavelength.
    /**
     * \returns the current value of m_lam_wfs.
     */
    realT lam_wfs();
 
    /// Set the number of pixels in the WFS for each WFS mode
    /** This sets the vector.
     */
    void npix_wfs( const std::vector<realT> &npix /**< [in] is the new values of m_npix_wfs */ );
 
    /// Set the number of pixels in the WFS for a given WFS mode
    /** This sets the entry in the vector, if it has the proper length.
     */
    void npix_wfs( int idx,   ///< [in] the index of the WFS mode
                   realT npix ///< [in] is the new value of m_npix_wfs[idx]
    );
 
    /// Get the number of pixels in the WFS for a given WFS mode
    /**
     * \returns the current value of m_npix_wfs[idx]
     */
    realT npix_wfs( size_t idx /**< [in] the index of the WFS mode.*/ );
 
    /// Get the number of pixels in the WFS for all WFS modes
    /**
     * \returns the current value of m_npix_wfs
     */
    std::vector<realT> npix_wfs();
 
    /// Set the value of the WFS readout noise for each WFS mode
    /** This sets the vector.
     */
    void ron_wfs( const std::vector<realT> &nron /**< [in] is the new value of m_ron_wfs */ );
 
    /// Set the value of the WFS readout noise for a given bining mode
    /** This sets the entry in the vector if the vector has the proper length
     */
    void ron_wfs( int idx,   ///< [in] the index of the WFS mode
                  realT nron ///< [in] is the new value of m_ron_wfs [idx]
    );
 
    /// Get the value of the WFS readout noise for a given WFS mode
    /**
     * \returns the current value of m_ron_wfs[idx]
     */
    realT ron_wfs( size_t idx /**< [in] the index of the WFS mode.*/ );
 
    /// Get the value of the WFS readout noise for all WFS modes
    /**
     * \returns the current value of m_ron_wfs
     */
    std::vector<realT> ron_wfs();
 
    /// Set the value of the background fluxes.
    /** This sets the value of the background flux for each WFS mode
     */
    void Fbg( const std::vector<realT> &fbg /**< [in] is the new value of m_Fbg */ );
 
    /// Set a single value of the background flux for a given WFS mode
    /** This sets the entry in the vector if the vector has the proper length
     */
    void Fbg( int idx,  ///< [in] the index of the WFS mode
              realT fbg ///< [in] is the new value of m_Fbg[idx]
    );
 
    /// Get the value of the background flux for a given WFS mode
    /**
     * \returns m_Fbg[idx]
     */
    realT Fbg( size_t idx /**< [in] the index of the WFS mode.*/ );
 
    /// Get the value of the background flux for all WFS modes
    /**
     * \returns the current value of m_Fbg
     */
    std::vector<realT> Fbg();
 
    /// Set the value of the minimum WFS exposure times.
    /** This sets the vector of the minimum WFS exposure times
     */
    void minTauWFS( const std::vector<realT> &ntau /**< [in] is the new value of m_minTauWFS */ );
 
    /// Set a single value of the minimum WFS exposure time for a given WFS mode.
    /** This sets the entry in the vector if the vector has sufficient length.
     */
    void minTauWFS( size_t idx, ///< [in] the index of the WFS mode
                    realT ntau  ///< [in] is the new value of m_minTauWFS
    );
 
    /// Get the value of the minimum WFS exposure time for a given binning.
    /**
     * \returns the current value of m_minTauWFS[idx].
     */
    realT minTauWFS( size_t idx /**< [in] the index of the WFS mode.*/ );
 
    /// Get the values of the minimum WFS exposure time.
    /**
     * \returns the current value of m_minTauWFS.
     */
    std::vector<realT> minTauWFS();
 
    /// Set the value of the pixel binning flag
    /**
     */
    void bin_npix( bool bnp /**< [in] is the new value of m_bin_npix */ );
 
    /// Get the value of the pixel binngin flag
    /**
     * \returns
     */
    bool bin_npix();
 
    /// Get the value of the optimum binning factor/index.
    /** If WFS modes are used, this is the mode index (0 to N-1).  If not, it is 1 minus the pixel binning factor.
     * It is always 1 minus the actuator binning factor.
     * This calls opt_d() to perform the optimization if needed.
     *
     * \returns the current value of m_bin_opt.
     */
    int bin_opt();
 
    /// Set the value of the WFS exposure time.
    /**
     */
    void tauWFS( realT ntau /**< [in] is the new value of m_tauWFS */ );
 
    /// Get the value of the minimum WFS exposure time.
    /**
     * \returns the current value of m_tauWFS.
     */
    realT tauWFS();
 
    /// Set the value of m_deltaTau.
    /**
     */
    void deltaTau( realT ndel /**< [in] is the new value of m_deltaTau*/ );
 
    /// Get the value of m_deltaTau.
    /**
     * \returns the current value of m_deltaTau.
     */
    realT deltaTau();
 
    /// Set the value of m_optTau.
    /**
     */
    void optTau( bool ot /**< [in] is the new value of m_optTau */ );
 
    /// Get the value of m_optTau.
    /**
     * \returns the current value of m_optTau.
     */
    bool optTau();
 
    /// Set the science wavelength.
    /**
     */
    void lam_sci( realT nlam /**< [in] is the new value of m_lam_sci */ );
 
    /// Get the science wavelength.
    /**
     * \returns the current value of m_lam_sci.
     */
    realT lam_sci();
 
    /// Set the zenith angle, and its secant.
    /**
     */
    void zeta( realT nz /**< [in] The new value of m_zeta */ );
 
    /// Get the zenith angle
    /**
     * \return the current value of m_zeta
     */
    realT zeta();
 
    /// Get the zecant of the zenith angle
    /**
     * \return the current value of m_secZeta
     */
    realT secZeta();
 
    /// Set the value of m_fit_mn_max
    /**
     */
    void fit_mn_max( int mnm /**< [in] is the new value of m_fit_mn_max */ );
 
    /// Get the value of m_fit_mn_max
    /**
     */
    int fit_mn_max();
 
    /// Set the value of the circularLimit flag
    /** If this is true, then the spatial frequency limits are treated as a circle
     * rather than a square.  This will create a circular dark hole.
     */
    void circularLimit( bool cl /**< [in] the new value of the circularLimit flag*/ );
 
    /// Get the value of the circularLimit flag
    /**
     * \returns the current value of m_circularLimit
     */
    bool circularLimit();
 
    /// Set the value of spatialFilter_ku
    /**
     */
    void spatialFilter_ku( realT ku /**< [in] is the new value of spatialFilter_ku*/ );
 
    /// Get the value of spatialFilter_ku
    /** \returns the current value of m_spatialFilter_ku
     */
    realT spatialFilter_ku();
 
    /// Set the value of spatialFilter_kv
    /**
     */
    void spatialFilter_kv( realT kv /**< [in] is the new value of spatialFilter_kv*/ );
 
    /// Get the value of spatialFilter_kv
    /** \returns the current value of m_spatialFilter_kv
     */
    realT spatialFilter_kv();
 
    /// Set the value of the non-common path WFE.
    /**
     */
    void ncp_wfe( realT nwfe /**< [in] is the new value of m_ncp_wfe*/ );
 
    /// Get the value of the non-common path WFE.
    /**
     * \returns the current value of m_ncp_wfe.
     */
    realT ncp_wfe();
 
    /// Set the value of the non-common path WFE PSD index.
    /**
     */
    void ncp_alpha( realT alpha /**< [in] is the new value of m_ncp_alpha*/ );
 
    /// Get the value of the non-common path WFE PSD index.
    /**
     * \returns the current value of m_ncp_alpha.
     */
    realT ncp_alpha();
 
    /// Set the value of the 0 magnitude photon rate
    /** This is the photon rate at the WFS, photons/sec.
     *
     */
    void F0( realT nF0 /**< [in] is the new value of m_F0.*/ );
 
    /// Get the value of the 0 magnitude photon rate
    /**
     * \returns the current value of m_F0.
     */
    realT F0();
 
    /// Set the value of the Star's magnitude
    /**
     */
    void starMag( realT nmag /**< [in] is the new value of m_starMag.*/ );
 
    /// Get the value of the Star's magnitude
    /**
     * \returns the current value of m_starMag
     */
    realT starMag();
 
    /// The photon flux at a given star magnitude.
    /**
     * \returns the photon flux of the star in the bandpass assumed by F0.
     */
    realT Fg( realT mag /**< [in] is the magnitude of the star. */ );
 
    /// Get the photon rate at the current Star magnitude.
    /** Calculates \f$ F_\gamma = F_0 10^{-0.4 m} \f$ where \f$ F_0 \f$ is m_F0 and \f$ m \f$ is m_starMag.
     *
     * \returns the current value of the current photon rate.
     */
    realT Fg();
 
    /// Access the array of controlledModes
    /** This array indicates which modes are controlled after optimizeStrehl chooses them.
     *
     * \returns a const reference to m_controlledModes.
     */
    const Eigen::Array<int, -1, -1> &controlledModes();
 
    /** \name Measurement Error
     * Calculating the WFE due to WFS measurement noise.
     * @{
     */
 
    /// Calculate the terms of the signal to noise ratio squared (S/N)^2 for the WFS measurement
    /** The S/N squared is
       \f[
       (S/N)^2 = \frac{ F_\gamma^2 \tau_{wfs}^2 }{ F_\gamma \tau_{wfs} + n_{pix} F_{bg} \tau_{wfs} + n_{pix}
      \sigma_{ron}^2 } \f]
 
      *
      * \returns the S/N squared
      */
    realT
    signal2Noise2( realT &Nph,    ///< [out] the number of photons
                   realT tau_wfs, ///< [in] specifies the WFS exposure time.  If 0, then optimumTauWFS is used
                   realT d,       ///< [in] the actuator spacing in meters, used if binning WFS pixels
                   int b ///< [in] the binning parameter.  Either the WFS mode index, or the binning factor minus 1.
    );
 
    /// Calculate the measurement noise at a spatial frequency and specified actuator spacing
    /** Calculates the wavefront phase variance due measurement noise at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
     *
     * \returns the measurement noise in rad^2 rms at the science wavelength
     */
    realT
    measurementError( realT m, ///< [in] specifies the u component of the spatial frequency
                      realT n, ///< [in] specifies the v component of the spatial frequency
                      realT d, ///< [in] the actuator spacing in meters
                      int b ///< [in] the binning parameter.  Either the WFS mode index, or the binning factor minus 1.
    );
 
    /// Calculate the total measurement error over all corrected spatial frequencies
    /** This totals the measurement WFE at the optimum actuator spacing and binning.
     *
     * \overload
     *
     * \returns the total WFE due to measurement error.
     */
    realT measurementErrorTotal();
 
    ///@}
 
    /** \name Time Delay Error
     * Calculating the WFE due to time delay.
     * @{
     */
 
    /// Calculate the time delay at a spatial frequency at the optimum exposure time and the specified actuator spacing
    /** Calculates the wavefront phase variance due to time delay at \f$ f = (m/D)\hat{u} + (n/D)\hat{v} \f$.
     *
     * \returns the measurement noise in rad^2 rms at the science wavelength
     */
    realT
    timeDelayError( realT m, ///< [in] specifies the u component of the spatial frequency
                    realT n, ///< [in] specifies the v component of the spatial frequency
                    realT d, ///< [in] the actuator spacing, in meters
                    int b    ///< [in] the binning parameter.  Either the WFS mode index, or the binning factor minus 1.
    );
 
    /// Calculate the time delay error over all corrected spatial frequencies
    /** This totals the time delay WFE at the optimum actuator spacing and binning.
     *
     * \overload
     *
     * \returns the total WFE due to time delay.
     */
    realT timeDelayErrorTotal();
 
    ///@}
 
    /** \name Fitting Error
     * Calculating the WFE due to uncorrected spatial frequencies.
     * @{
     */
 
    /// Calculate the fitting error at a specific spatial frequency.
    /**
     * \returns the fitting error in rad^2 rms at the science wavelength at (m,n).
     */
    realT fittingError( realT m, ///< [in] specifies the u component of the spatial frequency
                        realT n  ///< [in] specifies the v component of the spatial frequency
    );
 
    /// Calculate the total fitting error over all uncorrected spatial frequencies.
    /** This totals the fitting WFE for the optimum actuator spacing and binning.
     *
     * \overload
     *
     * \returns the total fitting error.
     */
    realT fittingErrorTotal();
 
    ///@}
 
    /** \name Chromatic Errors
     * Calculating the WFE due to chromaticity in scintillaion, index of refraction, and dispersion.
     * @{
     */
 
    /// Calculate the wavefront error due to scintillation chromaticity in the OPD at a spatial frequency.
    /** This totals the WFE from scintillation chromaticity in the OPD.
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromScintOPDError( int m, ///< [in] specifies the u component of the spatial frequency
                              int n  ///< [in] specifies the v component of the spatial frequency
    );
 
    /// Calculate the wavefront error due to scintillation chromaticity in the OPD over all spatial frequencies.
    /** This totals the WFE from scintillation chromaticity in the OPD for the optimum actuator spacing.
     *
     * \overload
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromScintOPDErrorTotal();
 
    /// Calculate the wavefront error due to scintillation chromaticity in amplitude for a given spatial frequency.
    /**
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromScintAmpError( int m, ///< [in] specifies the u component of the spatial frequency
                              int n  ///< [in] specifies the v component of the spatial frequency
    );
 
    /// Calculate the wavefront error due to scintillation chromaticity in amplitude over all spatial frequencies.
    /**
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromScintAmpError();
 
    /// Calculate the wavefront error due to chromaticity in the index of refraction at a given spatial frequency.
    /** This totals the WFE from chromaticity in the index of refraction.
     *
     * \overload
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromIndexError( int m, ///< [in] specifies the u component of the spatial frequency
                           int n  ///< [in] specifies the v component of the spatial frequency
    );
 
    /// Calculate the wavefront error due to chromaticity in the index of refraction at a specific spatial frequency.
    /** This totals the WFE from chromaticity in the index of refraction for the optimum actuator spacing.
     *
     * \overload
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT chromIndexErrorTotal();
 
    /// Calculate the wavefront error due to dispersive anisoplanatism in the OPD at a given spatial frequency
    /** This calculates the WFE from dispersive anisoplanatism in the OPD.
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT dispAnisoOPDError( int m, ///< [in] specifies the u component of the spatial frequency
                             int n  ///< [in] specifies the v component of the spatial frequency
    );
 
    /// Calculate the wavefront error due to dispersive anisoplanatism in the OPD over all specific spatial frequencies.
    /** This totals the WFE from dispersive anisoplanatism in the OPD for the optimum actuator spacing.
     *
     * \overload
     *
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT dispAnisoOPDErrorTotal();
 
    /// Calculate the wavefront error due to dispersive anisoplanatism in the amplitude at a specific spatial frequency.
    /**
     * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
     */
    realT dispAnisoAmpError();
 
    ///@}
 
    /** \name Optimum Parameters
     * Functions to calculate the optimum integration time and actuator spacing.
     *
     * @{
     */
 
    /// Calculate the optimum exposure time for a given spatial frequency at a specified actuator spacing
    /** Finds the optimum exposure time at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
     *
     * \todo Check inclusion of X in parameters
     *
     * \returns the optimum expsoure time.
     */
    realT
    optimumTauWFS( realT m, ///< [in] is the spatial frequency index in u
                   realT n, ///< [in] is the spatial frequency index in v
                   realT d, ///< [in] the actuator spacing, in meters
                   int b    ///< [in] the binning parameter.  Either the WFS mode index, or the binning factor minus 1.
    );
 
    /// Calculate the optimum exposure time for a given spatial frequency at the optimum actuator spacing.
    /** Finds the optimum exposure time at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
     * If optd == false this uses the minimum actuator spacing.
     *
     * \todo Check inclusion of X in parameters
     *
     * \returns the optimum expsoure time.
     */
    realT optimumTauWFS( realT m, ///< [in] is the spatial frequency index in u
                         realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the optimum actuator spacing.
    /** Finds the value of m_d_opt where the fitting error is less than than the combined time delay and measurement
     * error.
     *
     * \returns the current value of m_d_opt.
     */
    realT d_opt();
 
    /// @}
 
    /// Calculate the NCP variance at a spatial frequency.
    /** Finds the NCP variance at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
     *
     * \returns the NCP variance at k.
     */
    realT ncpError( int m, ///< [in] is the spatial frequency index in u
                    int n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the total NCP variance in rad^2.
    /**
     * \returns the total NCP variance.
     */
    realT ncpError();
 
    /** \name Overall WFE and Strehl Ratio
     * Functions to calculate the total WFE and Strehl ratio.
     *
     * @{
     */
 
  protected:
    /// Calculate the component WFE, total WFE, and Strehl ratio.
    /** This should only be called when something changes.
     */
    void calcStrehl();
 
  public:
    /// Get the current value of the total WFE variance.
    /** If no changes, merely returns m_wfeVar.  Calls calcStrehl if there are changes.
     *
     * \returns the current value of m_wfeVar;
     */
    realT wfeVar();
 
    /// Get the Strehl ratio for a given actuator pitch and WFS WFS mode.
    /**
      * Strehl is calculated using the extended Marechal approximation:
      *
      \f[
       S = e^{-\sigma_{wfe}^2}
       \f]
      * where \f$ \sigma_{wfe}^2 \f$ is the WFE for the spacing and binning
      *
      * \returns the strehl for these values
      */
    realT strehl( realT d, ///< [in] the actuator spacing, in meters
                  int b    ///< [in] the binning parameter.  Either the WFS mode index, or the binning factor minus 1.
    );
 
    /// Get the current Strehl ratio.
    /** If no changes, merely returns m_strehl.  Calls calcStrehl if there are changes.
      * Strehl is calculated using the extended Marechal approximation:
      *
      \f[
       S = e^{-\sigma_{wfe}^2}
       \f]
      * where \f$ \sigma_{wfe}^2 \f$ is the current value of m_wfeVar.
      *
      * \returns the current value of m_strehl.
      */
    realT strehl();
 
    /// @}
 
    /** \name Contrast
     * Calculating contrast
     *
     * @{
     */
 
    /// Worker function for raw contrast fuctions
    /** The logic in this calculation is the same regardless of component, except for the calculation of variance.
     * The varFunc function pointer is used to calculate the variance, otherwise  This handles the other details such
     * as Strehl normalization, fitting error (if appropriate), and bounds checking.
     *
     * \note this gives the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam varFuncT  is a function-pointer-to-member (of this class) with signature realT(*f)(realT, realT)
     */
    template <typename varFuncT>
    realT C_( int m,             ///< [in] is the spatial frequency index in u/
              int n,             ///< [in] is the spatial frequency index in v.
              bool normStrehl,   ///< [in] flag controls whether the contrast is normalized by Strehl ratio/
              varFuncT varFunc,  ///< [in] the variance function to use.
              int doFittingError ///< [in] flag to describe how fitting error should be considered for this term:
                                 ///< FITTING_ERROR_NO, FITTING_ERROR_ZERO, FITTING_ERROR_X, or FITTING_ERROR_Y.
    );
 
    /// Worker function for the contrast-map functions.
    /** The map calculation is the same for all terms, except for the calculation of variance at each pixel.
     * The Cfunc function pointer is used to actually get the contrast.
     *
     * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image
     * \tparam CfuncT is a function-pointer-to-member (of this class) with signature realT (*f)(realT, realT, bool)
     */
    template <typename imageT, typename CfuncT>
    void C_Map( imageT &im,     ///< [out] the map image to be filled in.
                CfuncT Cfunc,   ///< [in] the raw contrast function to use for filling in the map.
                bool normStrehl ///< [in] flag controlling whether or not the map is normalized by Strehl
    );
 
    /// Calculate the residual variance due to uncorrected phase at a spatial frequency.
    /** Used to calculate contrast \ref C0().
     *
     * \returns variance at (m,n).
     */
    realT C0var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to uncorrected phase, C0.
    /** Contrast C0 is the uncorrected phase, with the effects of scintillation included.  See Guyon (2005) \cite
     * guyon_2005, and the updated derivation in Males \& Guyon (2018) \cite males_guyon_2018.
     *
     * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \returns C0.
     */
    realT C0( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in] is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio
    );
 
    /// Calculate a 2D map of contrast C0
    /**
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C0Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to uncorrected amplitude at a spatial frequency.
    /** Used to calculate contrast \ref C1().
     *
     * \returns variance at (m,n).
     */
    realT C1var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to uncorrected amplitude, C1.
    /** Contrast C1 is the uncorrected amplitude due to scintillation.  See Guyon (2005) \cite guyon_2005, and the
     * updated derivation in Males \& Guyon (2018) \cite males_guyon_2018.
     *
     * \returns C0.
     */
    realT C1( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in] is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio
    );
 
    /// Calculate a 2D map of contrast C1.
    /** The contrast is Strehl-normalized here.
     *
     * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C1Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to measurement and time delay errors in phase/OPD at a spatial frequency.
    /** Used to calculate contrast \ref C2().
     *
     * \returns variance at (m,n).
     */
    realT C2var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to measurement and time delay errors in phase/OPD at a spatial frequency.
    /** Contrast C2 is just the total variance due to time delay and measurement errors,
     * divided by the Strehl ratio.   See Guyon (2005) \cite guyon_2005, and the updated
     * derivation in Males \& Guyon (2018) \cite males_guyon_2018.
     *
     * \returns C2.
     */
    realT C2( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in]  is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast \ref C2().
    /** The contrast is Strehl-normalized here.
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C2Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to measurement and time delay errors in amplitude at a spatial frequency.
    /** Used to calculate contrast \ref C3().
     *
     * \returns variance at (m,n).
     */
    realT C3var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to measurement and time delay errors in amplitude at a spatial frequency.
    /** Contrast C3 is just the total variance due to time delay and measurement errors,
     * divided by the Strehl ratio.    See Guyon (2005) \cite guyon_2005, and the updated
     * derivation in Males \& Guyon (2018) \cite males_guyon_2018.
     *
     * \returns C3.
     */
    realT C3( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in]  is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast C3.
    /** The contrast is Strehl-normalized here.
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C3Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to scintilation-OPD chromaticity.
    /** Used to calculate contrast \ref C4().
     *
     * \returns variance at (m,n).
     */
    realT C4var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to scintilation-OPD chromaticity.
    /** Contrast C4 is due to the chromaticity of scintillation, causing the ODP measurement to be slightly incorrect at
     * the science wavelength. See Guyon (2005) \cite guyon_2005, and the updated derivation in Males \& Guyon (2018)
     * \cite males_guyon_2018.
     *
     * \returns C4.
     */
    realT C4( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in]  is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast C4
    /** The contrast is Strehl-normalized here.
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C4Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to to scintilation-amplitude chromaticity.
    /** Used to calculate contrast \ref C5().
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \returns variance at (m,n).
     */
    realT C5var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to scintilation-amplitude chromaticity.
    /** Contrast C5 is due to the chromaticity of scintillation, causing the amplitude measurement to be slightly
     * incorrect at the science wavelength. See Guyon (2005) \cite guyon_2005, and the updated derivation in Males \&
     * Guyon (2018) \cite males_guyon_2018.
     *
     * \returns C4.
     */
    realT C5( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in]  is the spatial frequency index in v
              bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast C5
    /** The contrast is Strehl-normalized here.
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C5Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to chromaticity of the index of refraction of air.
    /** Used to calculate contrast \ref C6().
     *
     * \returns variance at (m,n).
     */
    realT C6var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to chromaticity of the index of refraction of air.
    /** Contrast C6 is due to the index of refraction of air being wavelength dependent,  causing the ODP measurement to
     * be slightly incorrect at the science wavelength.   See Guyon (2005) \cite guyon_2005, and the updated derivation
     * in Males \& Guyon (2018) \cite males_guyon_2018.
     *
     * \returns C6.
     */
    realT C6( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in] is the spatial frequency index in v
              bool normStrehl = true ///< flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast C6
    /** The contrast is Strehl-normalized here.
     *
     * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C6Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    /// Calculate the residual variance due to dispersive anisoplanatism.
    /** Used to calculate contrast \ref C7().
     *
     * \returns variance at (m,n).
     */
    realT C7var( realT m, ///< [in] is the spatial frequency index in u
                 realT n  ///< [in] is the spatial frequency index in v
    );
 
    /// Calculate the contrast due to dispersive anisoplanatism.
    /** Contrast C7 is due to atmospheric dispersion causing light at different wavelengths to take different paths
     * through atmospheric turbulence. See Fitzgerald (2017, in prep) \cite fitzgerald_2017
     *
     * \returns C7.
     */
    realT C7( realT m,               ///< [in] is the spatial frequency index in u
              realT n,               ///< [in] is the spatial frequency index in v
              bool normStrehl = true ///< flag controls whether the contrast is normalized by Strehl ratio.
    );
 
    /// Calculate a 2D map of contrast C7
    /** The contrast is Strehl-normalized here.
     *
     *  \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
     *
     * \tparam imageT is an Eigen-like image type.
     */
    template <typename imageT>
    void C7Map( imageT &map, ///< [in] the map image to be filled in with contrast
                bool normStrehl );
 
    ///@}
 
    /// Output current parameters to a stream
    /** Outputs a formatted list of all current parameters.
     *
     */
    iosT &dumpAOSystem( iosT &ios /**< [in] a std::ostream-like stream. */ );
 
    /// Setup the configurator to configure this class
    /**
     * todo: "\test 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
    /**
     * \todo: "\test 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, class inputSpectT, typename iosT>
aoSystem<realT, inputSpectT, iosT>::aoSystem()
{
    wfsBeta<wfs<realT, iosT>>( nullptr ); // allocate a default ideal WFS.
}
 
template <typename realT, class inputSpectT, typename iosT>
aoSystem<realT, inputSpectT, iosT>::~aoSystem()
{
    if( m_wfsBeta && m_ownWfsBeta )
    {
        delete m_wfsBeta;
    }
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::loadGuyon2005()
{
    atm.loadGuyon2005();
 
    F0( 1.75e9 * 0.25 * math::pi<realT>() * 64. * 0.18 ); // Converting to photons/sec
    lam_wfs( 0.55e-6 );
    lam_sci( 1.6e-6 );
    D( 8. );
    starMag( 5 );
 
    // The rest of Guyon 05 is very ideal
    npix_wfs( std::vector<realT>( { 12868 } ) );
    ron_wfs( std::vector<realT>( { 0.0 } ) );
    Fbg( std::vector<realT>( { 0.0 } ) );
 
    d_min( 8.0 / 1e3 );               // Allow super fine sampling
    minTauWFS( (realT)( 1. / 1e9 ) ); // Just be super fast.
    tauWFS( 1. / 1e9 );               // Just be super fast
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::loadMagAOX()
{
    atm.loadLCO();
    F0( 4.2e10 );
    lam_wfs( 0.791e-6 );
    lam_sci( 0.656e-6 );
 
    npix_wfs( std::vector<realT>( { 9024 } ) );
    ron_wfs( std::vector<realT>( { 0.57 } ) );
    Fbg( std::vector<realT>( { 0.22 } ) );
 
    d_min( std::vector<realT>( { 6.5 / 48.0 } ) );
    minTauWFS( std::vector<realT>( { 1. / 3622. } ) );
    tauWFS( 1. / 3622. );
 
    D( 6.5 );
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::loadGMagAOX()
{
 
    loadMagAOX();
 
    F0( 7.6e10 * 368.0 / ( 0.25 * math::pi<realT>() * 6.5 * 6.5 * ( 1. - 0.29 * 0.29 ) ) ); // Scaled up.
 
    D( 25.4 );
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::D( realT nD )
{
    m_D = nD;
    psd.D( m_D );
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::D()
{
    return m_D;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::d_min( const std::vector<realT> &nd )
{
    m_d_min.assign( nd.begin(), nd.end() );
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::d_min( int idx, realT nd )
{
    if( m_d_min.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_min", "idx larger than m_d_min" );
    }
 
    m_d_min[idx] = nd;
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::d_min( size_t idx )
{
    if( m_d_min.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_min", "idx larger than m_d_min" );
    }
 
    return m_d_min[idx];
}
 
template <typename realT, class inputSpectT, typename iosT>
std::vector<realT> aoSystem<realT, inputSpectT, iosT>::d_min()
{
    return m_d_min;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::optd( bool od )
{
    m_optd = od;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
bool aoSystem<realT, inputSpectT, iosT>::optd()
{
    return m_optd;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::optd_delta( realT odd )
{
    m_optd_delta = odd;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::optd_delta()
{
    return m_optd_delta;
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename wfsT>
void aoSystem<realT, inputSpectT, iosT>::wfsBeta( const wfsT &w )
{
    if( m_wfsBeta && m_ownWfsBeta )
    {
        delete m_wfsBeta;
        m_wfsBeta = nullptr;
    }
 
    m_wfsBeta = (wfs<realT, iosT> *)&w;
    m_ownWfsBeta = false;
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename wfsT>
void aoSystem<realT, inputSpectT, iosT>::wfsBeta( const wfsT *w )
{
    if( m_wfsBeta && m_ownWfsBeta )
    {
        delete m_wfsBeta;
        m_wfsBeta = nullptr;
    }
 
    if( w )
    {
        m_wfsBeta = (wfs<realT, iosT> *)w;
        m_ownWfsBeta = false;
    }
    else
    {
        m_wfsBeta = new wfsT;
        m_ownWfsBeta = true;
    }
}
 
template <typename realT, class inputSpectT, typename iosT>
wfs<realT, iosT> *aoSystem<realT, inputSpectT, iosT>::wfsBeta()
{
    return m_wfsBeta;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::beta_p( realT m, realT n )
{
    if( m_wfsBeta == 0 )
        mxThrowException( err::paramnotset, "aoSystem::beta_p", "The WFS is not assigned." );
 
    return m_wfsBeta->beta_p( m, n, m_D, d_opt(), atm.r_0( m_lam_sci ) );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::beta_r( realT m, realT n )
{
    if( m_wfsBeta == 0 )
        mxThrowException( err::paramnotset, "aoSystem::beta_r", "The WFS is not assigned." );
 
    return m_wfsBeta->beta_r( m, n, m_D, d_opt(), atm.r_0( m_lam_sci ) );
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::opticalGain( realT og )
{
    m_opticalGain = og;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::opticalGain()
{
    return m_opticalGain;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::lam_wfs( realT nlam )
{
    m_lam_wfs = nlam;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::lam_wfs()
{
    return m_lam_wfs;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::npix_wfs( const std::vector<realT> &npix )
{
    m_npix_wfs.resize( npix.size() );
    for( size_t n = 0; n < npix.size(); ++n )
    {
        m_npix_wfs[n] = npix[n];
    }
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::npix_wfs( int idx, realT npix )
{
    if( m_npix_wfs.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::npix_wfs", "idx larger than m_npix_wfs" );
    }
 
    m_npix_wfs[idx] = npix;
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::npix_wfs( size_t idx )
{
    if( m_npix_wfs.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::npix_wfs", "idx larger than m_npix_wfs" );
    }
 
    return m_npix_wfs[idx];
}
 
template <typename realT, class inputSpectT, typename iosT>
std::vector<realT> aoSystem<realT, inputSpectT, iosT>::npix_wfs()
{
    return m_npix_wfs;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::ron_wfs( const std::vector<realT> &nron )
{
    m_ron_wfs.resize( nron.size() );
    for( size_t n = 0; n < nron.size(); ++n )
    {
        m_ron_wfs[n] = nron[n];
    }
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::ron_wfs( int idx, realT nron )
{
    if( m_ron_wfs.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::ron_wfs", "idx larger than m_ron_wfs" );
    }
 
    m_ron_wfs[idx] = nron;
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::ron_wfs( size_t idx )
{
    if( m_ron_wfs.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::ron_wfs", "idx larger than m_ron_wfs" );
    }
 
    return m_ron_wfs[idx];
}
 
template <typename realT, class inputSpectT, typename iosT>
std::vector<realT> aoSystem<realT, inputSpectT, iosT>::ron_wfs()
{
    return m_ron_wfs;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::Fbg( const std::vector<realT> &fbg )
{
    m_Fbg.resize( fbg.size() );
    for( size_t n = 0; n < fbg.size(); ++n )
    {
        m_Fbg[n] = fbg[n];
    }
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::Fbg( int idx, realT fbg )
{
    if( m_Fbg.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSyste,::Fbg", "idx larger than m_Fbg" );
    }
 
    m_Fbg[idx] = fbg;
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::Fbg( size_t idx )
{
    if( m_Fbg.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSyste,::Fbg", "idx larger than m_Fbg" );
    }
 
    return m_Fbg[idx];
}
 
template <typename realT, class inputSpectT, typename iosT>
std::vector<realT> aoSystem<realT, inputSpectT, iosT>::Fbg()
{
    return m_Fbg;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::minTauWFS( const std::vector<realT> &ntau )
{
    m_minTauWFS.resize( ntau.size() );
    for( size_t n = 0; n < ntau.size(); ++n )
    {
        m_minTauWFS[n] = ntau[n];
    }
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::minTauWFS( size_t idx, realT ntau )
{
    if( m_minTauWFS.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::minTauWFS", "idx larger than m_ntau_wfs" );
    }
 
    m_minTauWFS[idx] = ntau;
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::minTauWFS( size_t idx )
{
    if( m_minTauWFS.size() < idx + 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::minTauWFS", "idx larger than m_ntau_wfs" );
    }
 
    return m_minTauWFS[idx];
}
 
template <typename realT, class inputSpectT, typename iosT>
std::vector<realT> aoSystem<realT, inputSpectT, iosT>::minTauWFS()
{
    return m_minTauWFS;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::bin_npix( bool bnp )
{
    if( bnp != m_bin_npix )
    {
        m_bin_npix = bnp;
        m_specsChanged = true;
        m_dminChanged = true;
    }
}
 
template <typename realT, class inputSpectT, typename iosT>
bool aoSystem<realT, inputSpectT, iosT>::bin_npix()
{
    return m_bin_npix;
}
 
template <typename realT, class inputSpectT, typename iosT>
int aoSystem<realT, inputSpectT, iosT>::bin_opt()
{
    d_opt();
    return m_bin_opt;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::tauWFS( realT ntau )
{
    m_tauWFS = ntau;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::tauWFS()
{
    return m_tauWFS;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::deltaTau( realT ndel )
{
    m_deltaTau = ndel;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::deltaTau()
{
    return m_deltaTau;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::optTau( bool ot )
{
    m_optTau = ot;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
bool aoSystem<realT, inputSpectT, iosT>::optTau()
{
    return m_optTau;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::lam_sci( realT nlam )
{
    m_lam_sci = nlam;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::lam_sci()
{
    return m_lam_sci;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::zeta( realT nz )
{
    m_zeta = nz;
    m_secZeta = 1 / cos( m_zeta );
 
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::zeta()
{
    return m_zeta;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::secZeta()
{
    return m_secZeta;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::fit_mn_max( int mnm )
{
    if( mnm < 0 )
        mnm = 0;
    m_fit_mn_max = mnm;
}
 
template <typename realT, class inputSpectT, typename iosT>
int aoSystem<realT, inputSpectT, iosT>::fit_mn_max()
{
    return m_fit_mn_max;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::circularLimit( bool cl )
{
    m_circularLimit = cl;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
bool aoSystem<realT, inputSpectT, iosT>::circularLimit()
{
    return m_circularLimit;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::spatialFilter_ku( realT kx )
{
    m_spatialFilter_ku = fabs( kx );
    m_specsChanged = true; // not sure if needed
    m_dminChanged = true;  // not sure if needed
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::spatialFilter_ku()
{
    return m_spatialFilter_ku;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::spatialFilter_kv( realT ky )
{
    m_spatialFilter_kv = fabs( ky );
    m_specsChanged = true; // not sure if needed
    m_dminChanged = true;  // not sure if needed
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::spatialFilter_kv()
{
    return m_spatialFilter_kv;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::ncp_wfe( realT nwfe )
{
    m_ncp_wfe = nwfe;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::ncp_wfe()
{
    return m_ncp_wfe;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::ncp_alpha( realT alpha )
{
    m_ncp_alpha = alpha;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::ncp_alpha()
{
    return m_ncp_alpha;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::F0( realT nF0 )
{
    m_F0 = nF0;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::F0()
{
    return m_F0;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::starMag( realT nmag )
{
    m_starMag = nmag;
    m_specsChanged = true;
    m_dminChanged = true;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::starMag()
{
    return m_starMag;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::Fg( realT mag )
{
    return m_F0 * pow( 10.0, -0.4 * mag );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::Fg()
{
    return Fg( m_starMag );
}
 
template <typename realT, class inputSpectT, typename iosT>
const Eigen::Array<int, -1, -1> &aoSystem<realT, inputSpectT, iosT>::controlledModes()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
 
    return m_controlledModes;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::signal2Noise2( realT &Nph, realT tau_wfs, realT d, int b )
{
    Nph = Fg() * tau_wfs;
 
    double binfact = 1.0;
    int binidx = 0;
    if( m_bin_npix )
    {
        if( m_npix_wfs.size() == 1 ) // Only if "true binning", otherwise the WFS mode vectors handle it.
        {
            binfact = 1. / pow( (realT)b + 1, 2 );
        }
        else
            binidx = b;
    }
 
    return pow( Nph, 2 ) / ( m_npix_wfs[binidx] * binfact * ( m_Fbg[binidx] * tau_wfs + pow( m_ron_wfs[binidx], 2 ) ) );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::measurementError( realT m, realT n, realT d, int b )
{
    if( m == 0 and n == 0 )
        return 0;
 
    realT tau_wfs;
 
    if( m_optTau )
        tau_wfs = optimumTauWFS( m, n, d, b );
    else
        tau_wfs = m_tauWFS;
 
    if( m_wfsBeta == 0 )
        mxThrowException( err::paramnotset, "aoSystem::measurementError", "The WFS is not assigned." );
 
    realT beta_p = m_wfsBeta->beta_p( m, n, m_D, d, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
    realT beta_r = m_wfsBeta->beta_r( m, n, m_D, d, atm.r_0( m_lam_wfs ) ); // sqrt(m_opticalGain);
    realT Nph = 0;
    realT snr2 = signal2Noise2( Nph, tau_wfs, d, b );
 
    return ( pow( beta_r, 2 ) / snr2 + pow( beta_p, 2 ) / Nph ) * pow( m_lam_wfs / m_lam_sci, 2 );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::measurementErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeMeasurement;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::timeDelayError( realT m, realT n, realT d, int b )
{
    if( m == 0 and n == 0 )
        return 0;
 
    realT k = sqrt( m * m + n * n ) / m_D;
 
    realT tau_wfs;
 
    if( m_optTau )
        tau_wfs = optimumTauWFS( m, n, d, b );
    else
        tau_wfs = m_tauWFS;
 
    realT tau = tau_wfs + m_deltaTau;
 
    ///\todo handle multiple layers
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           sqrt( atm.X( k, m_lam_sci, m_secZeta ) ) * pow( math::two_pi<realT>() * atm.v_wind() * k, 2 ) *
           pow( tau, 2 );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::timeDelayErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeTimeDelay;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::fittingError( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
 
    ///\todo handle multiple layers
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::fittingErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeFitting;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::chromScintOPDError( int m, int n )
{
    return C4var( m, n );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::chromScintOPDErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeChromScintOPD;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::chromScintAmpError( int m, int n )
{
    return C5var( m, n );
 
#if 0
   int mn_max = floor(0.5*m_D/d_opt());
   
   realT sum = 0;
 
   for(int m = -mn_max; m <= mn_max; ++m)
   {
      for(int n = -mn_max; n <= mn_max; ++n)
      {
         if(n == 0 && m == 0) continue;
         
         sum += C5var(m,n);
      }
   }
 
   return sum;
#endif
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::chromIndexError( int m, int n )
{
    return C6var( m, n );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::chromIndexErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeChromIndex;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::dispAnisoOPDError( int m, int n )
{
    return C7var( m, n );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::dispAnisoOPDErrorTotal()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
    return m_wfeAnisoOPD;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::dispAnisoAmpError()
{
    return 0;
 
#if 0
   int mn_max = floor(0.5*m_D/d_opt());
   
   realT sum = 0;
 
   for(int m = -mn_max; m <= mn_max; ++m)
   {
      for(int n = -mn_max; n <= mn_max; ++n)
      {
         if(n == 0 && m == 0) continue;
         
         if( m_circularLimit )
         {
            if( m*m + n*n > mn_max*mn_max ) continue;
         }
         
         sum += C8var(m,n);
      }
   }
 
   return sum;
#endif
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::optimumTauWFS(
    realT m,
    realT n,
    realT dact, // here called dact due to parameter collision with root-finding
    int bbin )
{
    if( m_D == 0 )
    {
        mxError( "aoSystem::optimumTauWFS", MXE_PARAMNOTSET, "Diameter (D) not set." );
        return -1;
    }
 
    if( m_F0 == 0 )
    {
        mxError( "aoSystem::optimumTauWFS", MXE_PARAMNOTSET, "0-mag photon flux (F0) not set." );
        return -1;
    }
 
    double binfact = 1.0;
    int binidx = 0; // index into WFS configurations
    if( m_bin_npix )
    {
        if( m_npix_wfs.size() > 1 )
        {
            binidx = bbin;
        }
        else
        {
            binfact = 1.0 / pow( (realT)( bbin + 1 ), 2 );
        }
    }
 
    realT k = sqrt( m * m + n * n ) / m_D;
 
    realT F = Fg();
 
    if( m_wfsBeta == 0 )
        mxThrowException( err::paramnotset, "aoSystem::beta_p", "The WFS is not assigned." );
 
    realT beta_p = m_wfsBeta->beta_p( m, n, m_D, dact, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
    realT beta_r = m_wfsBeta->beta_r( m, n, m_D, dact, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
 
    // Set up for root finding:
    realT a, b, c, d, e;
 
    ///\todo handle multiple layers
    realT Atmp = 2 * pow( atm.lam_0(), 2 ) * psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) *
                 ( atm.X( k, m_lam_wfs, m_secZeta ) ) * pow( math::two_pi<realT>() * atm.v_wind() * k, 2 );
 
    a = Atmp;
    b = Atmp * m_deltaTau;
    c = 0;
    d = -1 * ( pow( m_lam_wfs, 2 ) / F ) *
        ( ( m_npix_wfs[binidx] * binfact * m_Fbg[binidx] / F ) * pow( beta_r, 2 ) + pow( beta_p, 2 ) );
    e = -2 * pow( m_lam_wfs, 2 ) * ( m_npix_wfs[binidx] * binfact ) * pow( m_ron_wfs[binidx], 2 ) * pow( beta_r, 2 ) /
        pow( F, 2 );
 
    std::vector<std::complex<realT>> x;
 
    // Get the roots
    math::quarticRoots( x, a, b, c, d, e );
 
    // Now pick the largest positive real root
    realT tauopt = 0.0;
 
    for( int i = 0; i < 4; i++ )
    {
        if( real( x[i] ) > 0 && imag( x[i] ) == 0 && real( x[i] ) > tauopt )
            tauopt = real( x[i] );
    }
 
    if( tauopt < m_minTauWFS[binidx] )
        tauopt = m_minTauWFS[binidx];
 
    return tauopt;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::optimumTauWFS( realT m, realT n )
{
    d_opt(); // also gets m_bin_opt;
    return optimumTauWFS( m, n, m_d_opt, m_bin_opt );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::d_opt()
{
    if( !m_dminChanged )
        return m_d_opt;
 
    // Validate the modes
    ///\todo this should be in a separate valideModes function
    if( m_npix_wfs.size() < 1 )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_opt", "npix_wfs must have at least one entry" );
    }
 
    if( m_d_min.size() != m_npix_wfs.size() )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_opt", "d_min must be the same size as npix_wfs" );
    }
 
    if( m_ron_wfs.size() != m_npix_wfs.size() )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_opt", "ron_wfs must be the same size as npix_wfs" );
    }
 
    if( m_Fbg.size() != m_npix_wfs.size() )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_opt", "F_bg must be the same size as npix_wfs" );
    }
 
    if( m_minTauWFS.size() != m_npix_wfs.size() )
    {
        mxThrowException( err::sizeerr, "aoSystem::d_opt", "minTauWFS must be the same size as npix_wfs" );
    }
 
    ///\todo investigate why we tolerate m_d_min = 0 here.  Is this just a config error?
    if( m_d_min.size() == 1 && m_d_min[0] == 0 )
    {
        m_d_opt = 1e-50;
        m_bin_opt = 0;
        m_dminChanged = false;
        return m_d_opt;
    }
 
    if( !m_optd )
    {
        m_d_opt = m_d_min[0];
        m_bin_opt = 0;
        m_dminChanged = false;
 
        return m_d_opt;
    }
 
    realT best_d = 0;
 
    if( m_bin_npix )
    {
        if( m_npix_wfs.size() > 1 ) // Optimize over the WFS modes
        {
            int best_idx = 0;
            best_d = m_d_min[0];
 
            realT bestStrehlOverall = 0;
 
            for( int b = 0; b < m_npix_wfs.size(); ++b )
            {
                realT d = m_d_min[b];
                realT s = strehl( d, b );
                realT bestStrehl = s;
 
                // Find the actuator pitch at which AO error is less than fitting error at Nyquist
                while( d < m_D / 2 )
                {
                    d += m_d_min[b] / m_optd_delta;
                    s = strehl( d, b );
 
                    if( s < bestStrehl ) // Strehl got worse.  Can't be <= otherwise small m_optd_deltas will
                                         // immediately out
                    {
                        d -= m_d_min[b] / m_optd_delta; // go back to previous d
                        break;
                    }
                    else
                        bestStrehl = s;
                }
 
                // Check if this is optimum so far
                if( bestStrehl >= bestStrehlOverall ) // >= to favor fewer actuators
                {
                    bestStrehlOverall = bestStrehl;
                    best_idx = b;
                    best_d = d;
                }
            }
 
            m_bin_opt = best_idx;
        }
        else
        {
            realT d = m_d_min[0];
            m_bin_opt = 0;
 
            realT s = strehl( d, m_bin_opt );
            realT bestStrehl = s;
 
            while( d < m_D / 2 )
            {
                d += m_d_min[0] / m_optd_delta;
                m_bin_opt = d / m_d_min[0] - 1;
                if( m_bin_opt < 0 )
                    m_bin_opt = 0;
 
                s = strehl( d, m_bin_opt );
 
                if( s < bestStrehl ) // Strehl got worse.  Can't be <= otherwise small m_optd_deltas will immediately
                                     // out
                {
                    d -= m_d_min[0] / m_optd_delta; // go back to previous d
                    m_bin_opt = d / m_d_min[0] - 1;
                    if( m_bin_opt < 0 )
                        m_bin_opt = 0;
                    break;
                }
                else
                    bestStrehl = s;
            }
            best_d = d;
        }
    }
    else
    {
        realT d = m_d_min[0];
        m_bin_opt = 0;
 
        realT s = strehl( d, m_bin_opt );
        realT bestStrehl = s;
 
        // Find the actuator pitch which minimizes total error (AO error is less than fitting error at Nyquist)
        while( d < m_D / 2 )
        {
            d += m_d_min[0] / m_optd_delta;
 
            s = strehl( d, m_bin_opt );
            if( s < bestStrehl ) // Strehl got worse.  Can't be <= otherwise small m_optd_deltas will immediately out
            {
                d -= m_d_min[0] / m_optd_delta; // go back to previous d
                break;
            }
            else
                bestStrehl = s;
        }
        best_d = d;
    }
 
    m_d_opt = best_d;
    m_dminChanged = false;
 
    return m_d_opt;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::ncpError( int m, int n )
{
    if( m == 0 and n == 0 )
        return 0;
 
    realT k = sqrt( m * m + n * n ) / m_D;
 
    realT kmax = m_fit_mn_max / m_D;
    realT kmin = 1. / m_D;
 
    realT beta;
    if( m_ncp_alpha != 2 )
    {
        beta = ( m_ncp_alpha - 2 ) / ( math::two_pi<realT>() ) * 1. /
               ( pow( kmin, -m_ncp_alpha + 2 ) - pow( kmax, -m_ncp_alpha + 2 ) );
    }
    else
    {
        beta = 1. / ( math::two_pi<realT>() * log( kmax / kmin ) );
    }
    return beta * ncpError() * pow( k, -m_ncp_alpha ) * pow( kmin, 2 );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::ncpError()
{
    return pow( m_ncp_wfe, 2 ) * pow( math::two_pi<realT>() / m_lam_sci, 2 );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::strehl( realT d, int b )
{
 
    int mn_max = floor( 0.5 * m_D / d );
 
    realT wfeVar = 0;
 
    bool controlled;
    for( int m = -m_fit_mn_max; m <= m_fit_mn_max; ++m )
    {
        for( int n = -m_fit_mn_max; n <= m_fit_mn_max; ++n )
        {
            controlled = false;
            if( m_circularLimit )
            {
                if( m * m + n * n <= mn_max * mn_max )
                    controlled = true;
            }
            else
            {
                if( abs( m ) <= mn_max && abs( n ) <= mn_max )
                    controlled = true;
            }
 
            if( controlled )
            {
                realT wfeMeasurement = measurementError( m, n, d, b );
                realT wfeTimeDelay = timeDelayError( m, n, d, b );
                realT wfeChromScintOPD = chromScintOPDError( m, n );
                realT wfeChromIndex = chromIndexError( m, n );
                realT wfeAnisoOPD = dispAnisoOPDError( m, n );
 
                realT wfeFitting = fittingError( m, n );
 
                if( wfeFitting < wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD )
                {
                    wfeVar += wfeFitting;
                }
                else
                {
                    wfeVar += wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD;
                }
            }
            else
            {
                wfeVar += fittingError( m, n );
            }
        }
    }
 
    wfeVar += ncpError();
    return exp( -1 * wfeVar );
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::calcStrehl()
{
    realT d = d_opt(); // have to run this to get m_bin_opt set.
    int b = m_bin_opt;
 
    int mn_max = floor( 0.5 * m_D / d );
 
    m_wfeMeasurement = 0;
    m_wfeTimeDelay = 0;
    m_wfeChromScintOPD = 0;
    m_wfeChromIndex = 0;
    m_wfeAnisoOPD = 0;
 
    m_wfeFitting = 0;
 
    m_wfeNCP = 0;
 
    m_controlledModes.resize( 2 * m_fit_mn_max + 1, 2 * m_fit_mn_max + 1 );
    m_controlledModes.setZero();
 
    bool controlled;
    for( int m = -m_fit_mn_max; m <= m_fit_mn_max; ++m )
    {
        for( int n = -m_fit_mn_max; n <= m_fit_mn_max; ++n )
        {
            if( m == 0 && n == 0 )
            {
                continue;
            }
 
            controlled = false;
            if( m_circularLimit )
            {
                if( m * m + n * n <= mn_max * mn_max )
                    controlled = true;
            }
            else
            {
                if( abs( m ) <= mn_max && abs( n ) <= mn_max )
                    controlled = true;
            }
 
            if( controlled )
            {
                realT wfeMeasurement = measurementError( m, n, d, b );
                realT wfeTimeDelay = timeDelayError( m, n, d, b );
                realT wfeChromScintOPD = chromScintOPDError( m, n );
                realT wfeChromIndex = chromIndexError( m, n );
                realT wfeAnisoOPD = dispAnisoOPDError( m, n );
 
                realT wfeFitting = fittingError( m, n );
 
                if( wfeFitting < wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD )
                {
                    m_wfeFitting += wfeFitting;
                }
                else // it's worth controlling this mode.
                {
                    m_wfeMeasurement += wfeMeasurement;
                    m_wfeTimeDelay += wfeTimeDelay;
                    m_wfeChromScintOPD += wfeChromScintOPD;
                    m_wfeChromIndex += wfeChromIndex;
                    m_wfeAnisoOPD += wfeAnisoOPD;
                    m_controlledModes( m_fit_mn_max + m, m_fit_mn_max + n ) = 1;
                }
            }
            else
            {
                m_wfeFitting += fittingError( m, n );
            }
        }
    }
 
    m_wfeNCP = ncpError();
 
    m_wfeVar = m_wfeMeasurement + m_wfeTimeDelay + m_wfeFitting + m_wfeChromScintOPD + m_wfeChromIndex + m_wfeAnisoOPD +
               m_wfeNCP;
 
    m_strehl = exp( -1 * m_wfeVar );
 
    m_specsChanged = false;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::wfeVar()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
 
    return m_wfeVar;
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::strehl()
{
    if( m_specsChanged || m_dminChanged )
        calcStrehl();
 
    return m_strehl;
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename varFuncT>
realT aoSystem<realT, inputSpectT, iosT>::C_( int m, int n, bool normStrehl, varFuncT varFunc, int doFittingError )
{
    if( m == 0 && n == 0 )
        return 0;
 
    // we run this to optimize regardless of whether we use it
    //  if already optimized then this is a minor hit
    realT S = strehl();
 
    if( !normStrehl )
        S = 1;
 
    if( doFittingError != FITTING_ERROR_NO )
    {
        int mn_max = m_D / ( 2. * d_opt() );
 
        if( m_controlledModes( m_fit_mn_max + m, m_fit_mn_max + n ) == false )
        {
            if( doFittingError == FITTING_ERROR_ZERO )
                return 0;
 
            realT fe = fittingError( m, n );
 
            realT k = sqrt( m * m + n * n ) / m_D;
 
            if( doFittingError == FITTING_ERROR_X )
                fe *= ( atm.X( k, m_lam_sci, m_secZeta ) );
            else if( doFittingError == FITTING_ERROR_Y )
                fe *= ( atm.Y( k, m_lam_sci, m_secZeta ) );
            else
            {
                std::cerr << "Unknown doFittingError\n";
                exit( -1 );
            }
 
            return fe / S;
        }
    }
    // get if not doing fitting error or if inside control region:
 
    realT var = ( this->*varFunc )( m, n );
 
    return var / S;
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT, typename CfuncT>
void aoSystem<realT, inputSpectT, iosT>::C_Map( imageT &im, CfuncT Cfunc, bool normStrehl )
{
    int dim1 = im.rows();
    int dim2 = im.cols();
 
    int mc = 0.5 * ( dim1 - 1 );
    int nc = 0.5 * ( dim2 - 1 );
 
    for( int i = 0; i < dim1; ++i )
    {
        int m = i - mc;
 
        for( int j = 0; j < dim2; ++j )
        {
            int n = j - nc;
 
            im( i, j ) = ( this->*Cfunc )( m, n, normStrehl );
        }
    }
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C0var( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
    ///\todo handle multiple layers
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           ( atm.X( k, m_lam_sci, m_secZeta ) );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C0( realT m, realT n, bool normStrehl )
{
 
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C0var, FITTING_ERROR_NO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C0Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C0, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C1var( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
 
    ///\todo handle multiple layers
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           ( atm.Y( k, m_lam_sci, m_secZeta ) );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C1( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C1var, FITTING_ERROR_NO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C1Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C1, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C2var( realT m, realT n )
{
    d_opt();
    return measurementError( m, n, d_opt(), m_bin_opt ) + timeDelayError( m, n, d_opt(), m_bin_opt );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C2( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C2var, FITTING_ERROR_X );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C2Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C2, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C3var( realT m, realT n )
{
    return 0; // measurementError(m, n) + timeDelayError(m,n);
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C3( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C3var, FITTING_ERROR_ZERO ); // FITTING_ERROR_Y);
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C3Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C3, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C4var( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
 
    ///\todo handle multiple layers
    // This does not need to be divided by X b/c we haven't multiplied by it, this is is C0/X.
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           atm.dX( k, m_lam_sci, m_lam_wfs );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C4( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C4var, FITTING_ERROR_ZERO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C4Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C4, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C5var( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
 
    ///\todo handle multiple layers
    // This does not need to be divided by Y b/c we haven't multiplied by it, this is is C1/Y.
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           atm.dY( k, m_lam_sci, m_lam_wfs );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C5( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C5var, FITTING_ERROR_ZERO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C5Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C5, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C6var( realT m, realT n )
{
    realT ni = atm.n_air( m_lam_sci );
    realT nw = atm.n_air( m_lam_wfs );
 
    return C0var( m, n ) * pow( ( ni - nw ) / ( ni - 1 ), 2 ); // Fixes error in Eqn 29
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C6( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C6var, FITTING_ERROR_ZERO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C6Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C6, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C7var( realT m, realT n )
{
    realT k = sqrt( m * m + n * n ) / m_D;
 
    ///\todo this needs to handle multiple layers -- do we violate an assumption of L_0 isn't constant?
    return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
           atm.X_Z( k, m_lam_wfs, m_lam_sci, m_secZeta );
}
 
template <typename realT, class inputSpectT, typename iosT>
realT aoSystem<realT, inputSpectT, iosT>::C7( realT m, realT n, bool normStrehl )
{
    return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C7var, FITTING_ERROR_ZERO );
}
 
template <typename realT, class inputSpectT, typename iosT>
template <typename imageT>
void aoSystem<realT, inputSpectT, iosT>::C7Map( imageT &im, bool normStrehl )
{
    C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C7, normStrehl );
}
 
template <typename realT, class inputSpectT, typename iosT>
iosT &aoSystem<realT, inputSpectT, iosT>::dumpAOSystem( iosT &ios )
{
    ios << "# AO Params:\n";
    ios << "#    D = " << D() << '\n';
 
    if( m_d_min.size() > 0 )
    {
        ios << "#    d_min = " << d_min( (size_t)0 );
        for( size_t n = 1; n < m_d_min.size(); ++n )
            ios << ',' << d_min( n );
        ios << '\n';
    }
    else
        ios << "#    d_min = null\n";
 
    ios << "#    optd = " << std::boolalpha << m_optd << '\n';
    ios << "#    d_opt_delta = " << optd_delta() << '\n';
    ios << "#    lam_sci = " << lam_sci() << '\n';
    ios << "#    F0 = " << F0() << '\n';
    ios << "#    starMag = " << starMag() << '\n';
    ios << "#    lam_sci = " << lam_sci() << '\n';
    ios << "#    zeta    = " << zeta() << '\n';
    ios << "#    lam_wfs = " << lam_wfs() << '\n';
 
    if( npix_wfs().size() > 0 )
    {
        ios << "#    npix_wfs = " << npix_wfs( (size_t)0 );
        for( size_t n = 1; n < npix_wfs().size(); ++n )
            ios << ',' << npix_wfs( n );
        ios << '\n';
    }
    else
        ios << "#    npix_wfs = null\n";
 
    if( ron_wfs().size() > 0 )
    {
        ios << "#    ron_wfs = " << ron_wfs( (size_t)0 );
        for( size_t n = 1; n < ron_wfs().size(); ++n )
            ios << ',' << ron_wfs( n );
        ios << '\n';
    }
    else
        ios << "#    ron_wfs = null\n";
 
    if( Fbg().size() > 0 )
    {
        ios << "#    Fbg = " << Fbg( (size_t)0 );
        for( size_t n = 1; n < Fbg().size(); ++n )
            ios << ',' << Fbg( n );
        ios << '\n';
    }
    else
        ios << "#    Fbg = null\n";
 
    if( minTauWFS().size() > 0 )
    {
        ios << "#    minTauWFS = " << minTauWFS( (size_t)0 );
        for( size_t n = 1; n < minTauWFS().size(); ++n )
            ios << ',' << minTauWFS( n );
        ios << '\n';
    }
    else
        ios << "#    minTauWFS = null\n";
 
    ios << "#    bin_npix = " << std::boolalpha << m_bin_npix << '\n';
    ios << "#    tauWFS = " << tauWFS() << '\n';
    ios << "#    optTau = " << std::boolalpha << m_optTau << '\n';
    ios << "#    deltaTau = " << deltaTau() << '\n';
    ios << "#    fit_mn_max = " << m_fit_mn_max << '\n';
    ios << "#    spatialFilter_ku = " << m_spatialFilter_ku << '\n';
    ios << "#    spatialFilter_kv = " << m_spatialFilter_kv << '\n';
    ios << "#    ncp_wfe = " << m_ncp_wfe << '\n';
    ios << "#    ncp_alpha = " << m_ncp_alpha << '\n';
 
    if( m_wfsBeta == 0 )
        mxThrowException( err::paramnotset, "aoSystem::beta_p", "The WFS is not assigned." );
    m_wfsBeta->dumpWFS( ios );
    psd.dumpPSD( ios );
    atm.dumpAtmosphere( ios );
 
    dumpGitStatus( ios );
 
    return ios;
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::setupConfig( app::appConfigurator &config )
{
    using namespace mx::app;
 
    // AO System configuration
    config.add( "aosys.wfs",
                "",
                "aosys.wfs",
                argType::Required,
                "aosys",
                "wfs",
                false,
                "string",
                "The WFS type: idealWFS, unmodPyWFS, asympModPyWFS, shwfs, calculatedWFS" );
    config.add( "aosys.wfs_beta_p",
                "",
                "aosys.wfs_beta_p",
                argType::Required,
                "aosys",
                "wfs_beta_p",
                false,
                "string",
                "The beta_p file path for calcualtedWFS" );
    config.add( "aosys.wfs_beta_r",
                "",
                "aosys.wfs_beta_r",
                argType::Required,
                "aosys",
                "wfs_beta_r",
                false,
                "string",
                "The beta_r file path for calcualtedWFS" );
    config.add( "aosys.wfs_sensitivity",
                "",
                "aosys.wfs_sensitivity",
                argType::Required,
                "aosys",
                "wfs_sensitivity",
                false,
                "bool",
                "Flag indicating that beta_p/beta_r are sensitivities (inverse) [default false]" );
    config.add(
        "aosys.D", "", "aosys.D", argType::Required, "aosys", "D", false, "real", "The telescope diameter [m]" );
    config.add( "aosys.d_min",
                "",
                "aosys.d_min",
                argType::Required,
                "aosys",
                "d_min",
                false,
                "real",
                "The minimum actuator spacing [m]" );
    config.add( "aosys.optd",
                "",
                "aosys.optd",
                argType::Optional,
                "aosys",
                "optd",
                false,
                "bool",
                "Whether or not the actuator spacing is optimized" );
    config.add( "aosys.optd_delta",
                "",
                "aosys.optd_delta",
                argType::Required,
                "aosys",
                "optd_delta",
                false,
                "bool",
                "The fractional change from d_min used in optimization.  Set to 1 (default) for integer binnings, > 1 "
                "for finer sampling." );
    config.add( "aosys.F0",
                "",
                "aosys.F0",
                argType::Required,
                "aosys",
                "F0",
                false,
                "real",
                "Zero-mag photon flux, [photons/sec]" );
    config.add( "aosys.lam_wfs",
                "",
                "aosys.lam_wfs",
                argType::Required,
                "aosys",
                "lam_wfs",
                false,
                "real",
                "WFS wavelength [m]" );
    config.add( "aosys.npix_wfs",
                "",
                "aosys.npix_wfs",
                argType::Required,
                "aosys",
                "npix_wfs",
                false,
                "vector<real>",
                "The number of pixels in the WFS" );
    config.add( "aosys.ron_wfs",
                "",
                "aosys.ron_wfs",
                argType::Required,
                "aosys",
                "ron_wfs",
                false,
                "vector<real>",
                "WFS readout noise [photons/read]" );
    config.add( "aosys.bin_npix",
                "",
                "aosys.bin_npix",
                argType::Required,
                "aosys",
                "bin_npix",
                false,
                "bool",
                "Whether or not WFS pixels are re-binned along with actuator spacing optimization" );
    config.add( "aosys.Fbg",
                "",
                "aosys.Fbg",
                argType::Required,
                "aosys",
                "Fbg",
                false,
                "vector<real>",
                "Background counts, [counts/pix/sec]" );
    config.add( "aosys.tauWFS",
                "",
                "aosys.tauWFS",
                argType::Required,
                "aosys",
                "tauWFS",
                false,
                "real",
                "WFS integration time [s]" );
    config.add( "aosys.minTauWFS",
                "",
                "aosys.minTauWFS",
                argType::Required,
                "aosys",
                "minTauWFS",
                false,
                "vector<real>",
                "Minimum WFS integration time [s]" );
    config.add( "aosys.deltaTau",
                "",
                "aosys.deltaTau",
                argType::Required,
                "aosys",
                "deltaTau",
                false,
                "real",
                "Loop delay [s]" );
    config.add( "aosys.optTau",
                "",
                "aosys.optTau",
                argType::Optional,
                "aosys",
                "optTau",
                false,
                "bool",
                "Whether or not the integration time is optimized" );
    config.add( "aosys.lam_sci",
                "",
                "aosys.lam_sci",
                argType::Required,
                "aosys",
                "lam_sci",
                false,
                "real",
                "Science wavelength [m]" );
    config.add(
        "aosys.zeta", "", "aosys.zeta", argType::Required, "aosys", "zeta", false, "real", "Zenith distance [rad]" );
    config.add( "aosys.fit_mn_max",
                "",
                "aosys.fit_mn_max",
                argType::Required,
                "aosys",
                "fit_mn_max",
                false,
                "real",
                "Maximum spatial frequency index to use for analysis" );
    config.add( "aosys.circularLimit",
                "",
                "aosys.circularLimit",
                argType::Optional,
                "aosys",
                "circularLimit",
                false,
                "bool",
                " Flag to indicate that the spatial frequency limit is circular, not square." );
    config.add( "aosys.spatialFilter_ku",
                "",
                "aosys.spatialFilter_ku",
                argType::Required,
                "aosys",
                "spatialFilter_ku",
                false,
                "real",
                "Spatial filter cutoff frequency in u [m^-1]" );
    config.add( "aosys.spatialFilter_kv",
                "",
                "aosys.spatialFilter_kv",
                argType::Required,
                "aosys",
                "spatialFilter_kv",
                false,
                "real",
                "Spatial filter cutoff frequency in v [m^-1]" );
    config.add( "aosys.ncp_wfe",
                "",
                "aosys.ncp_wfe",
                argType::Required,
                "aosys",
                "ncp_wfe",
                false,
                "real",
                "NCP WFE between 1 lambda/D and fit_mn_max [rad^2]" );
    config.add( "aosys.ncp_alpha",
                "",
                "aosys.ncp_alpha",
                argType::Required,
                "aosys",
                "ncp_alpha",
                false,
                "real",
                "PSD index for NCP WFE" );
    config.add(
        "aosys.starMag", "", "aosys.starMag", argType::Required, "aosys", "starMag", false, "real", "Star magnitude" );
    config.add( "aosys.starMags",
                "",
                "aosys.starMags",
                argType::Required,
                "aosys",
                "starMags",
                false,
                "real vector",
                "A vector of star magnitudes" );
 
    atm.setupConfig( config );
    psd.setupConfig( config );
}
 
template <typename realT, class inputSpectT, typename iosT>
void aoSystem<realT, inputSpectT, iosT>::loadConfig( app::appConfigurator &config )
{
    // WFS
    if( config.isSet( "aosys.wfs" ) )
    {
        std::string wfsStr;
        config( wfsStr, "aosys.wfs" );
 
        if( wfsStr == "ideal" )
        {
            wfsBeta<wfs<realT, iosT>>( nullptr );
        }
        else if( wfsStr == "unmodPyWFS" )
        {
            wfsBeta<pywfsUnmod<realT, iosT>>( nullptr );
        }
        else if( wfsStr == "asympModPyWFS" )
        {
            wfsBeta<pywfsModAsymptotic<realT>>( nullptr );
        }
        else if( wfsStr == "SHWFS" )
        {
            wfsBeta<shwfs<realT, iosT>>( nullptr );
        }
        else if( wfsStr == "calculatedWFS" )
        {
            wfsBeta<calculatedWFS<realT, iosT>>( nullptr );
 
            calculatedWFS<realT, iosT> *cwfs = static_cast<calculatedWFS<realT, iosT> *>( m_wfsBeta );
            config( cwfs->m_beta_p_file, "aosys.wfs_beta_p" );
            config( cwfs->m_beta_r_file, "aosys.wfs_beta_r" );
            bool sens = cwfs->m_sensitivity;
            config( sens, "aosys.wfs_sensitivity" );
            if( config.isSet( "aosys.wfs_sensitivity" ) )
            {
                cwfs->m_sensitivity = sens;
            }
        }
        else
        {
            mxThrowException( err::invalidconfig, "aoSystem::loadConfig", "unknown WFS " + wfsStr + " specified" );
        }
    }
 
    // We load the default value, get the config value (which will be the default if not set), and
    // then if it was set, call the setter function so any side effects are captured.
 
    // diameter
    realT nD = D();
    config( nD, "aosys.D" );
    if( config.isSet( "aosys.D" ) )
        D( nD );
 
    // d_min
    std::vector<realT> nd_min = d_min();
    config( nd_min, "aosys.d_min" );
    if( config.isSet( "aosys.d_min" ) )
        d_min( nd_min );
 
    bool noptd = optd();
    config( noptd, "aosys.optd" );
    if( config.isSet( "aosys.optd" ) )
        optd( noptd );
 
    realT noptd_delta = optd_delta();
    config( noptd_delta, "aosys.optd_delta" );
    if( config.isSet( "aosys.optd_delta" ) )
        optd_delta( noptd_delta );
 
    realT nlam_wfs = lam_wfs();
    config( nlam_wfs, "aosys.lam_wfs" );
    if( config.isSet( "aosys.lam_wfs" ) )
        lam_wfs( nlam_wfs );
 
    // npix_wfs
    std::vector<realT> nnpix_wfs = npix_wfs();
    config( nnpix_wfs, "aosys.npix_wfs" );
    if( config.isSet( "aosys.npix_wfs" ) )
        npix_wfs( nnpix_wfs );
 
    // ron_wfs
    std::vector<realT> nron_wfs = ron_wfs();
    config( nron_wfs, "aosys.ron_wfs" );
    if( config.isSet( "aosys.ron_wfs" ) )
        ron_wfs( nron_wfs );
 
    // Fbg
    std::vector<realT> nFbg = Fbg();
    config( nFbg, "aosys.Fbg" );
    if( config.isSet( "aosys.Fbg" ) )
        Fbg( nFbg );
 
    // minTauWFS
    std::vector<realT> nminTauWFS = minTauWFS();
    config( nminTauWFS, "aosys.minTauWFS" );
    if( config.isSet( "aosys.minTauWFS" ) )
        minTauWFS( nminTauWFS );
 
    // bin_npix
    bool nbin_npix = bin_npix();
    config( nbin_npix, "aosys.bin_npix" );
    if( config.isSet( "aosys.bin_npix" ) )
        bin_npix( nbin_npix );
 
    // tauWFS
    realT ntauWFS = tauWFS();
    config( ntauWFS, "aosys.tauWFS" );
    if( config.isSet( "aosys.tauWFS" ) )
        tauWFS( ntauWFS );
 
    // deltaTau
    realT ndeltaTau = deltaTau();
    config( ndeltaTau, "aosys.deltaTau" );
    if( config.isSet( "aosys.deltaTau" ) )
        deltaTau( ndeltaTau );
 
    // optTau
    bool noptTau = optTau();
    config( noptTau, "aosys.optTau" );
    if( config.isSet( "aosys.optTau" ) )
        optTau( noptTau );
 
    // lam_sci
    realT nlam_sci = lam_sci();
    config( nlam_sci, "aosys.lam_sci" );
    if( config.isSet( "aosys.lam_sci" ) )
        lam_sci( nlam_sci );
 
    // zeta
    realT nzeta = zeta();
    config( nzeta, "aosys.zeta" );
    if( config.isSet( "aosys.zeta" ) )
        zeta( nzeta );
 
    // fit_mn_max
    realT fmnm = fit_mn_max();
    config( fmnm, "aosys.fit_mn_max" );
    if( config.isSet( "aosys.fit_mn_max" ) )
        fit_mn_max( fmnm );
 
    // circularLimit
    bool cl = circularLimit();
    config( cl, "aosys.circularLimit" );
    if( config.isSet( "aosys.circularLimit" ) )
        circularLimit( cl );
 
    // spatialFilter_ku
    realT ku = spatialFilter_ku();
    config( ku, "aosys.spatialFilter_ku" );
    if( config.isSet( "aosys.spatialFilter_ku" ) )
        spatialFilter_ku( ku );
 
    // spatialFilter_kv
    realT kv = spatialFilter_kv();
    config( kv, "aosys.spatialFilter_kv" );
    if( config.isSet( "aosys.spatialFilter_kv" ) )
        spatialFilter_kv( kv );
 
    // ncp_wfe
    realT nwfe = ncp_wfe();
    config( nwfe, "aosys.ncp_wfe" );
    if( config.isSet( "aosys.ncp_wfe" ) )
        ncp_wfe( nwfe );
 
    // ncp_alpha
    realT na = ncp_alpha();
    config( na, "aosys.ncp_alpha" );
    if( config.isSet( "aosys.ncp_alpha" ) )
        ncp_alpha( na );
 
    // F0
    realT nF0 = F0();
    config( nF0, "aosys.F0" );
    if( config.isSet( "aosys.F0" ) )
        F0( nF0 );
 
    // star_mag
    realT smag = starMag();
    config( smag, "aosys.starMag" );
    if( config.isSet( "aosys.starMag" ) )
        starMag( smag );
 
    atm.loadConfig( config );
    psd.loadConfig( config );
}
 
extern template class aoSystem<float, vonKarmanSpectrum<float>, std::ostream>;
 
extern template class aoSystem<double, vonKarmanSpectrum<double>, std::ostream>;
 
extern template class aoSystem<long double, vonKarmanSpectrum<long double>, std::ostream>;
 
#ifdef HASQUAD
extern template class aoSystem<__float128, vonKarmanSpectrum<__float128>, std::ostream>;
#endif
 
} // namespace analysis
} // namespace AO
} // namespace mx
 
#endif // aoSystem_hpp