From fd8be18dafb6b3cb1785b94dc4b86bbf7733edcb Mon Sep 17 00:00:00 2001 From: Jared Males Date: Sun, 9 Aug 2026 09:41:53 -0700 Subject: [PATCH] added theoretical zernike vars; added half_root_three --- bg/calcConstants.cpp | 9 +- include/math/constants.hpp | 96 ++++++++++++++- include/sigproc/zernike.hpp | 237 ++++++++++++++++++++++++++++++++++++ 3 files changed, 338 insertions(+), 4 deletions(-) diff --git a/bg/calcConstants.cpp b/bg/calcConstants.cpp index b61ccce55..c93ce3d98 100644 --- a/bg/calcConstants.cpp +++ b/bg/calcConstants.cpp @@ -22,6 +22,7 @@ int main() floatT twosqrt2log2 = static_cast( 2.0 ) * sqrt( static_cast( 2.0 ) * log( static_cast( 2.0 ) ) ); + std::cout << " "; // label buffer for( int i = 0; i < 10; ++i ) // print digit scale { @@ -31,7 +32,9 @@ int main() } std::cout << "1\n"; - std::cout << "tan_arcsec: " << tan_arcsec << "\n"; - std::cout << "twosqrt2log2: " << twosqrt2log2 << "\n"; - std::cout << "ln_two: " << ln_two() << "\n"; + std::cout << "tan_arcsec: " << tan_arcsec << "\n"; + std::cout << "twosqrt2log2: " << twosqrt2log2 << "\n"; + std::cout << "ln_two: " << ln_two() << "\n"; + std::cout << "root_three: " << sqrt(static_cast(3)) << '\n'; + std::cout << "half_root_three: " << sqrt(static_cast(3))/static_cast(2) << '\n'; } diff --git a/include/math/constants.hpp b/include/math/constants.hpp index 47f6905c1..ae98868cd 100644 --- a/include/math/constants.hpp +++ b/include/math/constants.hpp @@ -37,13 +37,20 @@ namespace mx { namespace math { -// 0 1 2 3 4 5 6 7 8 9 1 +// Constants to 100 digits for casting +// 1 +// 1 2 3 4 5 6 7 8 9 0 +// 01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 #define MX_INTERNAL_PI_100 \ ( 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 ) #define MX_INTERNAL_ROOT2_100 \ ( 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727 ) #define MX_INTERNAL_LN2_100 \ ( 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875 ) +#define MX_INTERNAL_ROOT3_100 \ + ( 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576 ) +#define MX_INTERNAL_HALF_ROOT3_100 \ + ( 0.8660254037844386467637231707529361834714026269051903140279034897259665084544000185405730933786242878 ) /// Get the value of pi /** Specializations provided for float, double, long double, and quad if supported. Can default to boost for other types @@ -316,6 +323,93 @@ constexpr __float128 ln_two<__float128>() } #endif +/// Get the value of sqrt(3) +/** Specializations provided for float, double, long double, and quad (if supported). Can default to boost for other + * types if MX_INCLUDE_BOOST is defined. + * + * \ingroup genconstants + */ +template +constexpr T root_three() +{ +#ifdef MX_INCLUDE_BOOST + return boost::math::constants::root_three(); +#else + static_assert( + std::is_fundamental::value || !std::is_fundamental::value, + "root_three not specialized for type T, and MX_INCLUDE_BOOST is not defined, so I can't just use boost." ); + return 0; +#endif +} + +template <> +constexpr float root_three() +{ + return static_cast( MX_INTERNAL_ROOT3_100 ); +} + +template <> +constexpr double root_three() +{ + return static_cast( MX_INTERNAL_ROOT3_100 ); +} + +template <> +constexpr long double root_three() +{ + return static_cast( MX_INTERNAL_ROOT3_100 ); +} + +#ifdef HASQUAD +template <> +constexpr __float128 root_three<__float128>() +{ + return static_cast<__float128>( MX_INTERNAL_ROOT3_100 ); +} +#endif + +/// Get the value of sqrt(3)/2 +/** Specializations provided for float, double, long double, and quad (if supported). + * + * Note there is no boost equivalent. + * + * \ingroup genconstants + */ +template +constexpr T half_root_three() +{ + static_assert( + std::is_fundamental::value || !std::is_fundamental::value, + "half_root_three not specialized for type T, and there is no boost equivalent." ); + return 0; +} + +template <> +constexpr float half_root_three() +{ + return static_cast( MX_INTERNAL_HALF_ROOT3_100 ); +} + +template <> +constexpr double half_root_three() +{ + return static_cast( MX_INTERNAL_HALF_ROOT3_100 ); +} + +template <> +constexpr long double half_root_three() +{ + return static_cast( MX_INTERNAL_HALF_ROOT3_100 ); +} + +#ifdef HASQUAD +template <> +constexpr __float128 half_root_three<__float128>() +{ + return static_cast<__float128>( MX_INTERNAL_HALF_ROOT3_100 ); +} +#endif + /// Get the value of 1/3 /** Wrapper for boost constant. Specializations provided for float, double, and long double. * diff --git a/include/sigproc/zernike.hpp b/include/sigproc/zernike.hpp index 6b8d49df4..671d76569 100644 --- a/include/sigproc/zernike.hpp +++ b/include/sigproc/zernike.hpp @@ -758,6 +758,243 @@ realT zernikePTrefoil( const realT &kD /**< [in] Spatial frequency in diameter u return 32 * pow( math::func::jincN( 4, math::pi() * kD ), 2 ); } +/// Get the degrees of correction coefficient for Zernike polynomials in Kolmogorov turbulence. +/** Returns the coefficient from Table IV of Noll (1976) \cite noll_1976 for the given index. + * Given this, the total variance in radians at wavelength $\lambda$ for a diameter $D$ + * aperture after correcting $j$ modes can be calculated from + * \f[ + * realT c = zernikeModeDOCKolmogorov( j ); + * var_lambda = c * pow(D/r_0_lambda, math::five_thirds) + * \f] + * \returns 0 if noll_j is 0, indicating an error + * \returns the coefficient + */ +template +realT zernikeModeDOCKolmogorov( unsigned noll_j /**< [in] the mode index, must be greater than 0 */ ) +{ + realT c; + + switch( noll_j ) + { + case 1: + c = 1.0299; + break; + case 2: + c = 0.582; + break; + case 3: + c = 0.134; + break; + case 4: + c = 0.111; + break; + case 5: + c = 0.0880; + break; + case 6: + c = 0.0648; + break; + case 7: + c = 0.0587; + break; + case 8: + c = 0.0525; + break; + case 9: + c = 0.0463; + break; + case 10: + c = 0.0401; + break; + case 11: + c = 0.0377; + break; + case 12: + c = 0.0352; + break; + case 13: + c = 0.0328; + break; + case 14: + c = 0.0304; + break; + case 15: + c = 0.0279; + break; + case 16: + c = 0.0267; + break; + case 17: + c = 0.0255; + break; + case 18: + c = 0.0243; + break; + case 19: + c = 0.0232; + break; + case 20: + c = 0.0220; + break; + case 21: + c = 0.0208; + break; + default: + if( noll_j == 0 ) + { + return 0; + } + + c = 0.2944 * pow( static_cast( noll_j ), -1 * math::half_root_three() ); + } + + return c; +} + +/// Get the degrees of correction for Zernike polynomials in Kolmogorov turbulence. +/** Returns the degree of correction from Table IV of Noll (1976) \cite noll_1976 for the given index. + * This is the total variance in radians for Fried parameter $r_0$ for a diameter $D$. + * Equivalent to: + * \f[ + * realT c = zernikeModeDOCKolmogorov( j ); + * var_lambda = c * pow(D/r_0_lambda, math::five_thirds) + * \f] + * \returns 0 if noll_j is 0, indicating an error + * \returns the coefficient + */ +template +realT zernikeModeDOCKolmogorov( unsigned noll_j, /**< [in] the mode index, must be greater than 0 */ + realT D, /**< [in] the aperture diameter, in same units as r_0 */ + realT r_0 /** <[in] the Fried parameter, in same units as D */ +) +{ + if( noll_j == 0 ) + { + return 0; + } + + return zernikeModeDOCKolmogorov( noll_j ) * pow( D / r_0, math::five_thirds() ); +} + +/// Get the difference in degrees of correction coefficient for Zernike polynomials in Kolmogorov turbulence. +/** Returns the difference in coefficients from Table IV of Noll (1976) \cite noll_1976 for the given index + * from the previous mode. Given this, the variance in radians-squared at wavelength $\lambda$ for a + * diameter $D$ aperture in the $j$-th mode can be calculated from + * \f[ + * realT c = zernikeModeDOCDiffKolmogorov( j ); + * var_j_lambda = c * pow(D/r_0_lambda, math::five_thirds) + * \f] + * + * \returns 0 if noll_j is 0 or 1, indicating an error + * \returns the coefficient difference + */ +template +realT zernikeModeDOCDiffKolmogorov( unsigned noll_j /**< [in] the mode index, must be greater than 1 */ ) +{ + realT c; + + switch( noll_j ) + { + case 2: + c = 1.0299 - 0.582; + break; + case 3: + c = 0.582 - 0.134; + break; + case 4: + c = 0.134 - 0.111; + break; + case 5: + c = 0.111 - 0.0880; + break; + case 6: + c = 0.0880 - 0.0648; + break; + case 7: + c = 0.0648 - 0.0587; + break; + case 8: + c = 0.0587 - 0.0525; + break; + case 9: + c = 0.0525 - 0.0463; + break; + case 10: + c = 0.0463 - 0.0401; + break; + case 11: + c = 0.0401 - 0.0377; + break; + case 12: + c = 0.0377 - 0.0352; + break; + case 13: + c = 0.0352 - 0.0328; + break; + case 14: + c = 0.0328 - 0.0304; + break; + case 15: + c = 0.0304 - 0.0279; + break; + case 16: + c = 0.0279 - 0.0267; + break; + case 17: + c = 0.0267 - 0.0255; + break; + case 18: + c = 0.0255 - 0.0243; + break; + case 19: + c = 0.0243 - 0.0232; + break; + case 20: + c = 0.0232 - 0.0220; + break; + case 21: + c = 0.0220 - 0.0208; + break; + case 22: + c = 0.0208 - 0.2944 * pow( static_cast( 22 ), -1 * math::half_root_three() ); + break; + default: + if( noll_j < 2 ) + { + return 0; + } + + c = 0.2944 * ( pow( static_cast( noll_j - 1 ), -1 * math::half_root_three() ) - + pow( static_cast( noll_j ), -1 * math::half_root_three() ) ); + } + + return c; +} + +/// Get the variance for a single Zernike polynomial in Kolmogorov turbulence. +/** Returns the variance from Table IV of Noll (1976) \cite noll_1976 for the given mode index. + * For the $j$-th mode on a diameter $D$ aperture and Fried parameter $r_0$ this is equivalent to calculating + * \f[ + * realT c = zernikeModeDOCDiffKolmogorov( j ); + * var_j_lambda = c * pow(D/r_0_lambda, math::five_thirds) + * \f] + * + * \returns 0 if noll_j is 0 or 1, indicating an error + * \returns the variance for the \param noll_j mode in radians squared. + */ +template +realT zernikeModeDOCDiffKolmogorov( unsigned noll_j, /**< [in] the mode index, must be greater than 0 */ + realT D, /**< [in] the aperture diameter, in same units as r_0 */ + realT r_0 /** <[in] the Fried parameter, in same units as D */ ) +{ + if( noll_j < 2 ) + { + return 0; + } + + return zernikeModeDOCDiffKolmogorov( noll_j ) * pow( D / r_0, math::five_thirds() ); +} + ///@} signal_processing } // namespace sigproc