Line data Source code
1 : /** \file aoPSDs.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Spatial power spectra used in adaptive optics.
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2016-2018 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef aoPSDs_hpp
28 : #define aoPSDs_hpp
29 :
30 : #include <string>
31 :
32 : #include "../../mxlib.hpp"
33 : #include "../../math/constants.hpp"
34 : #include "../../math/func/jinc.hpp"
35 :
36 : #include "aoConstants.hpp"
37 :
38 : #include "aoAtmosphere.hpp"
39 :
40 : namespace mx
41 : {
42 : namespace AO
43 : {
44 : namespace analysis
45 : {
46 :
47 : namespace PSDComponent
48 : {
49 : /// Enum to specify which component of the PSD to calcualte
50 : enum
51 : {
52 : phase, ///< The phase or OPD
53 : amplitude, ///< The amplitude
54 : dispPhase, ///< The phase component of dispersive anisoplanatism
55 : dispAmplitude ///< The amplitude component of dispersive anisoplanatism
56 : };
57 :
58 : std::string compName( int cc );
59 :
60 : int compNum( const std::string &name );
61 : } // namespace PSDComponent
62 :
63 : /// Manage calculations using the von Karman spatial power spectrum.
64 : /** This very general PSD has the form
65 : *
66 : * \f[
67 : \mathcal{P}(k) = \frac{\beta}{ ( k^2 + k_0^2) ^{\alpha/2}}
68 : \f]
69 : *
70 : * Where \f$ k \f$ is spatial frequency, \f$ \beta \f$ is a normalization constant, and outer scale \f$ L_0 \f$ is
71 : included via \f$ k_0 = 1/L_0 \f$.
72 : * For \f$ L_0 \rightarrow \infty \f$ this becomes a simple power law with index \f$ \alpha \f$.
73 : *
74 : * For atmospheric turbulence \f$ \alpha = 11/3 \f$ and the normalization is
75 : \f[
76 : \beta = \frac{\mathcal{A}_p}{r_0^{5/3}}
77 : \f]
78 : * where \f$ A_p = 0.0218...\f$ (see \ref mx::AO::constants::a_PSD ).
79 : *
80 : * \ingroup mxAOAnalytic
81 : */
82 : template <typename realT>
83 : struct vonKarmanSpectrum
84 : {
85 : typedef verbose::d verboseT;
86 :
87 : protected:
88 : bool m_subPiston{ true }; ///< flag controlling whether piston is subtracted from the PSD. Default is true.
89 : bool m_subTipTilt{
90 : false }; ///< flag controlling whether tip and tilt are subtracted from the PSD. Default is false.
91 :
92 : bool m_scintillation{ false }; ///< flag controlling whether or not scintillation is included
93 : int m_component{ PSDComponent::phase }; ///< If m_scintillation is true, this controls whether phase (0), amplitude
94 : ///< (1), or dispersive contrast (2) is returned.
95 :
96 : realT m_D{ 1.0 }; ///< Diameter used for piston and tip/tilt subtraction, in m. Default is 1 m.
97 :
98 : const char *m_id = "von Karman";
99 :
100 : // realT m_alpha {eleven_thirds<realT>()}; ///< The power-law index, 11/3 for Kolmogorov
101 :
102 : public:
103 : /// Default Constructor
104 : vonKarmanSpectrum();
105 :
106 : /// Constructor specifying the parameters
107 : /**
108 : */
109 : vonKarmanSpectrum( bool subP, ///< [in] is the value of m_subPiston.
110 : bool subT, ///< [in] is the value of m_subTipTilt.
111 : realT D ///< [in] is the value of m_D.
112 : );
113 :
114 : /// Get the value of m_subPiston
115 : /**
116 : * \returns m_subPiston
117 : */
118 : bool subPiston();
119 :
120 : /// Set the value of m_subPiston
121 : /**
122 : */
123 : void subPiston( bool sp /**< [in] is the new value of m_subPiston */ );
124 :
125 : /// Get the value of m_subTipTilt
126 : /**
127 : * \returns the value of m_subTipTilt.
128 : */
129 : bool subTipTilt();
130 :
131 : /// Set the value of the m_subTipTilt flag.
132 : /**
133 : */
134 : void subTipTilt( bool st /**< [in] the new value of m_subTipTilt.*/ );
135 :
136 : /// Get the value of m_scintillation
137 : /**
138 : * \returns m_scintillation
139 : */
140 : bool scintillation();
141 :
142 : /// Set the value of m_scintillation
143 : /**
144 : */
145 : void scintillation( bool sc /**< [in] the new value of m_scintillation*/ );
146 :
147 : /// Get the value of m_component
148 : /**
149 : * \returns m_component
150 : */
151 : int component();
152 :
153 : /// Set the value of m_component
154 : /**
155 : */
156 : int component( int cc /**< [in] the new value of m_component */ );
157 :
158 : /// Get the value of the diameter m_D.
159 : /**
160 : * \returns the current value of m_D, the diameter in m.
161 : */
162 : realT D();
163 :
164 : /// Set the aperture diameter
165 : /**
166 : */
167 : void D( realT nd /**< [in] the new diameter in m */ );
168 :
169 : /// Get the value of the PSD at spatial frequency k and a zenith distance.
170 : /**
171 : * \returns the von Karman PSD at the specified spatial frequency.
172 : * \returns -1 if an error occurs.
173 : *
174 : */
175 : template <class psdParamsT>
176 : realT operator()( psdParamsT &par, ///< [in] gives the PSD parameters.
177 : size_t layer_i,
178 : realT k, ///< [in] is the spatial frequency in m^-1.
179 : realT sec_zeta ///< [in] is the secant of the zenith distance.
180 : );
181 :
182 : /// Get the value of the PSD at spatial frequency k and wavelength lambda, and a zenith distance, with a WFS at a
183 : /// different wavelength
184 : /**
185 : *
186 : * \returns the von Karman PSD at the specified spatial frequency for the specified wavelength.
187 : * \returns -1 if an error occurs.
188 : */
189 : template <class psdParamsT>
190 : realT operator()( psdParamsT &par, ///< [in] gives the PSD parameters.
191 : size_t layer_i,
192 : realT k, ///< [in] is the spatial frequency in m^-1.
193 : realT lambda, ///< [in] is the observation wavelength in m
194 : realT lambda_wfs, ///< [in] is the wavefront measurement wavelength in m
195 : realT secZeta ///< [in] is the secant of the zenith distance
196 : );
197 :
198 : /// Get the fitting error for an actuator spacing d.
199 : /**
200 : * \todo generatlize for different alpha and beta
201 : *
202 : * \param atm gives the atmosphere parameters r_0 and L_0
203 : * \param d is the actuator spacing in m
204 : */
205 : // realT fittingError(aoAtmosphere<realT> &atm, realT d);
206 :
207 : template <typename iosT>
208 : iosT &dumpPSD( iosT &ios );
209 :
210 : /// Setup the configurator to configure this class
211 : void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
212 :
213 : /// Load the configuration of this class from a configurator
214 : void loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
215 : };
216 :
217 : template <typename realT>
218 52 : vonKarmanSpectrum<realT>::vonKarmanSpectrum()
219 : {
220 52 : }
221 :
222 : template <typename realT>
223 0 : vonKarmanSpectrum<realT>::vonKarmanSpectrum( bool subP, // [in] is the value of m_subPiston.
224 : bool subT, // [in] is the value of m_subTipTilt.
225 : realT D // [in] is the diameter
226 0 : )
227 : {
228 0 : m_subPiston = subP;
229 0 : m_subTipTilt = subT;
230 0 : m_D = D;
231 0 : }
232 :
233 : template <typename realT>
234 41076 : bool vonKarmanSpectrum<realT>::subPiston()
235 : {
236 41076 : return m_subPiston;
237 : }
238 :
239 : template <typename realT>
240 0 : void vonKarmanSpectrum<realT>::subPiston( bool sp /* [in] is the new value of m_subPiston */ )
241 : {
242 0 : m_subPiston = sp;
243 0 : }
244 :
245 : template <typename realT>
246 41076 : bool vonKarmanSpectrum<realT>::subTipTilt()
247 : {
248 41076 : return m_subTipTilt;
249 : }
250 :
251 : template <typename realT>
252 0 : void vonKarmanSpectrum<realT>::subTipTilt( bool st /* [in] the new value of m_subTipTilt */ )
253 : {
254 0 : m_subTipTilt = st;
255 0 : }
256 :
257 : template <typename realT>
258 4 : bool vonKarmanSpectrum<realT>::scintillation()
259 : {
260 4 : return m_scintillation;
261 : }
262 :
263 : template <typename realT>
264 0 : void vonKarmanSpectrum<realT>::scintillation( bool sc /* [in] the new value of m_scintillation */ )
265 : {
266 0 : m_scintillation = sc;
267 0 : }
268 :
269 : template <typename realT>
270 4 : int vonKarmanSpectrum<realT>::component()
271 : {
272 4 : return m_component;
273 : }
274 :
275 : template <typename realT>
276 14 : int vonKarmanSpectrum<realT>::component( int cc /* [in] the new value of m_component */ )
277 : {
278 14 : if( cc != PSDComponent::phase && cc != PSDComponent::amplitude && cc != PSDComponent::dispPhase &&
279 : cc != PSDComponent::dispAmplitude )
280 : {
281 0 : internal::mxlib_error_report<verboseT>( error_t::invalidarg, "Unknown component" );
282 : // mxError( "vonKarmanSpectrum::component", MXE_INVALIDARG, "Unknown component" );
283 0 : return -1;
284 : }
285 :
286 14 : m_component = cc;
287 14 : return 0;
288 : }
289 :
290 : template <typename realT>
291 4 : realT vonKarmanSpectrum<realT>::D()
292 : {
293 4 : return m_D;
294 : }
295 :
296 : template <typename realT>
297 39 : void vonKarmanSpectrum<realT>::D( realT nd /**< [in] the new diameter in m */ )
298 : {
299 39 : m_D = nd;
300 39 : }
301 :
302 : template <typename realT>
303 : template <class psdParamsT>
304 96450 : realT vonKarmanSpectrum<realT>::operator()( psdParamsT &par, //< [in] gives the PSD parameters.
305 : size_t layer_i,
306 : realT k, // [in] is the spatial frequency in m^-1.
307 : realT sec_zeta // [in] is the secant of the zenith distance.
308 : )
309 : {
310 : realT k02;
311 :
312 : ///\todo this needs to handle layers with different L_0
313 96450 : if( par.L_0( layer_i ) > 0 )
314 : {
315 96180 : k02 = ( 1 ) / ( par.L_0( layer_i ) * par.L_0( layer_i ) );
316 : }
317 : else
318 270 : k02 = 0;
319 :
320 96450 : if( k02 == 0 && k == 0 )
321 : {
322 0 : return 0;
323 : }
324 :
325 : realT Ppiston, Ptiptilt;
326 :
327 96450 : if( ( m_subPiston || m_subTipTilt ) )
328 : {
329 96450 : if( m_D == 0 )
330 : {
331 0 : internal::mxlib_error_report<verboseT>( error_t::paramnotset,
332 0 : "Diameter D not set for Piston and/or TT subtraction." );
333 :
334 : // mxError( "aoAtmosphere", MXE_PARAMNOTSET, "Diameter D not set for Piston and/or TT subtraction." );
335 0 : return -1;
336 : }
337 96450 : if( m_subPiston )
338 : {
339 96450 : Ppiston = pow( 2 * math::func::jinc( math::pi<realT>() * k * m_D ), 2 );
340 : }
341 : else
342 0 : Ppiston = 0;
343 :
344 96450 : if( m_subTipTilt )
345 : {
346 0 : Ptiptilt = pow( 4 * math::func::jincN( 2, math::pi<realT>() * k * m_D ), 2 );
347 : }
348 : else
349 96450 : Ptiptilt = 0;
350 : }
351 : else
352 : {
353 0 : Ppiston = 0;
354 0 : Ptiptilt = 0;
355 : }
356 :
357 96450 : return ( par.beta( layer_i ) * pow( k * k + k02, -1 * par.alpha( layer_i ) / 2 ) + par.beta_0( layer_i ) ) *
358 96450 : ( 1.0 - Ppiston - Ptiptilt ) * sec_zeta;
359 : }
360 :
361 : template <typename realT>
362 : template <class psdParamsT>
363 0 : realT vonKarmanSpectrum<realT>::operator()(
364 : psdParamsT &par, // [in] gives the PSD parameters.
365 : size_t layer_i,
366 : realT k, // [in] is the spatial frequency in m^-1.
367 : realT lambda, // [in] is the observation wavelength in m. Not used if par.nonKolmogorov==true
368 : realT lambda_wfs, // [in] is the wavefront measurement wavelength in m. Only used if m_scintillation==true.
369 : realT secZeta // [in] is the secant of the zenith distance.
370 : )
371 : {
372 0 : realT psd = operator()( par, layer_i, k, secZeta );
373 :
374 0 : if( par.nonKolmogorov() == false )
375 : {
376 0 : psd *= pow( par.lam_0() / lambda, 2 );
377 : }
378 :
379 0 : if( psd < 0 )
380 0 : return -1;
381 :
382 0 : if( m_scintillation )
383 : {
384 0 : if( m_component == PSDComponent::phase )
385 : {
386 0 : psd *= ( par.X( k, lambda, secZeta ) );
387 : }
388 0 : else if( m_component == PSDComponent::amplitude )
389 : {
390 0 : psd *= ( par.Y( k, lambda, secZeta ) );
391 : }
392 0 : else if( m_component == PSDComponent::dispPhase )
393 : {
394 0 : psd *= ( par.X_Z( k, lambda, lambda_wfs, secZeta ) );
395 : }
396 0 : else if( m_component == PSDComponent::dispAmplitude )
397 : {
398 0 : internal::mxlib_error_report<verboseT>( error_t::notimpl, "Dispersive-aniso amplitude not implemented" );
399 : // mxError( "vonKarmanSpectrum::operator()", MXE_NOTIMPL, "Dispersive-aniso amplitude not implemented" );
400 0 : return 0;
401 : }
402 : else
403 : {
404 0 : internal::mxlib_error_report<verboseT>( error_t::invalidarg, "Invalid component specified" );
405 : // mxError( "vonKarmanSpectrum::operator()", MXE_INVALIDARG, "Invalid component specified" );
406 0 : return 0;
407 : }
408 : }
409 :
410 0 : return psd;
411 : }
412 :
413 : /*template< typename realT>
414 : realT vonKarmanSpectrum<realT>::fittingError(aoAtmosphere<realT> &atm, realT d)
415 : {
416 : realT k0;
417 : if(atm.L_0(0) > 0)
418 : {
419 : k0 = 1/ (atm.L_0(0)*atm.L_0(0));
420 : }
421 : else k0 = 0;
422 :
423 : return (math::pi<realT>() * math::six_fifths<realT>())* constants::a_PSD<realT>()/ pow(atm.r_0(),
424 : math::five_thirds<realT>()) * (1./pow( pow(0.5/d,2) + k0, math::five_sixths<realT>()));
425 : }*/
426 :
427 : template <typename realT>
428 : template <typename iosT>
429 1 : iosT &vonKarmanSpectrum<realT>::dumpPSD( iosT &ios )
430 : {
431 1 : ios << "# PSD Parameters:" << '\n';
432 1 : ios << "# ID = " << m_id << '\n';
433 1 : ios << "# D = " << m_D << '\n';
434 1 : ios << "# subPiston = " << std::boolalpha << m_subPiston << '\n';
435 1 : ios << "# subTipTilt = " << std::boolalpha << m_subTipTilt << '\n';
436 1 : ios << "# Scintillation = " << std::boolalpha << m_scintillation << '\n';
437 1 : ios << "# Component = " << PSDComponent::compName( m_component ) << '\n';
438 1 : return ios;
439 : }
440 :
441 : template <typename realT>
442 14 : void vonKarmanSpectrum<realT>::setupConfig( app::appConfigurator &config )
443 : {
444 : using namespace mx::app;
445 :
446 14 : config.add( "psd.D",
447 28 : "",
448 28 : "psd.D",
449 : argType::Required,
450 28 : "psd",
451 28 : "D",
452 : false,
453 28 : "real",
454 42 : "Aperture diameter. Used for piston and tip/tilt subtraction." );
455 196 : config.add( "psd.subPiston", "", "psd.subPiston", argType::Required, "psd", "subPiston", false, "real", "" );
456 196 : config.add( "psd.subTipTilt", "", "psd.subTipTilt", argType::Required, "psd", "subTipTilt", false, "real", "" );
457 14 : config.add( "psd.scintillation",
458 28 : "",
459 28 : "psd.scintillation",
460 : argType::Required,
461 28 : "psd",
462 28 : "scintillation",
463 : false,
464 28 : "real",
465 42 : "" );
466 182 : config.add( "psd.component", "", "psd.component", argType::Required, "psd", "component", false, "real", "" );
467 14 : }
468 :
469 : template <typename realT>
470 14 : void vonKarmanSpectrum<realT>::loadConfig( app::appConfigurator &config )
471 : {
472 : // Here "has side effecs" means that the set function does more than simply copy the value.
473 :
474 28 : config( m_D, "psd.D" );
475 28 : config( m_subPiston, "psd.subPiston" );
476 28 : config( m_subTipTilt, "psd.subTipTilt" );
477 14 : config( m_scintillation, "psd.scintillation" );
478 :
479 14 : std::string cn = PSDComponent::compName( m_component );
480 14 : config( cn, "psd.component" );
481 14 : component( PSDComponent::compNum( cn ) );
482 14 : }
483 :
484 : extern template struct vonKarmanSpectrum<float>;
485 :
486 : extern template struct vonKarmanSpectrum<double>;
487 :
488 : extern template struct vonKarmanSpectrum<long double>;
489 :
490 : #ifdef HASQUAD
491 : extern template struct vonKarmanSpectrum<__float128>;
492 : #endif
493 :
494 : } // namespace analysis
495 : } // namespace AO
496 : } // namespace mx
497 :
498 : #endif // aoPSDs_hpp
|
