Line data Source code
1 : /** \file aoAtmosphere.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Provides a class to specify atmosphere parameters.
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : #ifndef aoAtmosphere_hpp
9 : #define aoAtmosphere_hpp
10 :
11 : #include <iostream>
12 : #include <numeric>
13 : #include <cmath>
14 : #include <cstdlib>
15 : #include <utility>
16 : #include <vector>
17 : #include <algorithm>
18 :
19 : #include "../../mxlib.hpp"
20 :
21 : #include "aoConstants.hpp"
22 :
23 : #include "../../math/constants.hpp"
24 : #include "../../math/floatUtils.hpp"
25 :
26 : #include "../../app/appConfigurator.hpp"
27 :
28 : namespace mx
29 : {
30 : namespace AO
31 : {
32 : namespace analysis
33 : {
34 :
35 : /// A class to specify atmosphere parameters and perform related calculations.
36 : /**
37 : * \todo layer outer scales need work
38 : * \todo layer PSD params isn't finished. Need to push that to PSD itself, and manage things like wavelength
39 : * dependence.
40 : *
41 : * \tparam realT is the real floating type in which all calculations are performed.
42 : *
43 : * \ingroup mxAOAnalytic
44 : */
45 : template <typename _realT>
46 : class aoAtmosphere
47 : {
48 : public:
49 : typedef _realT realT; ///< The real floating type in which all calculations are performed.
50 :
51 : /// Constructor
52 : aoAtmosphere();
53 :
54 : /// Validate the complete atmosphere configuration for calculation.
55 : /** Incremental setters may temporarily leave layer vectors incomplete. Call this at configuration and calculation
56 : * boundaries before using indexed or derived atmosphere values.
57 : *
58 : * \returns `error_t::noerror` when all scalar and layer invariants are satisfied, or a typed configuration error.
59 : */
60 : error_t validate() const;
61 :
62 : protected:
63 : realT m_r_0{ 0 }; ///< Fried's parameter, in m
64 :
65 : realT m_lam_0{ 0.5e-6 }; ///< Wavelength of Fried's parameter, in m
66 :
67 : std::vector<realT> m_layer_Cn2; ///< Vector of layer strengths.
68 :
69 : std::vector<realT> m_L_0; ///< The outer scale, in m
70 :
71 : std::vector<realT> m_l_0; ///< The inner scale of each layer, in m
72 :
73 : bool m_nonKolmogorov{ false }; ///< Flag indicating if non-Kolmogorov PSD parameters are used.
74 :
75 : std::vector<realT> m_beta{ 1 }; ///< The PSD normalization when in non-Kolmogorov mode.
76 :
77 : std::vector<realT> m_alpha{ 0 }; ///< The PSD exponent when in non-Kolmogorov mode.
78 :
79 : std::vector<realT> m_beta_0{ 0 }; ///< The PSD constant when in non-Kolmogorov mode.
80 :
81 : std::vector<realT> m_layer_z; ///< Vector of layer heights, in m, above the observatory.
82 :
83 : realT m_h_obs{ 0 }; ///< Height of the observatory above sea level, in m.
84 :
85 : realT m_H{ 8000 }; ///< The atmospheric scale height, in m.
86 :
87 : std::vector<realT> m_layer_v_wind; ///< Vector of layer wind speeds, in m/s.
88 :
89 : std::vector<realT> m_layer_dir; ///< Vector of layer wind directions, in radians.
90 :
91 : bool m_v_wind_updated{ false }; ///< whether or not m_v_wind has been updated after changes
92 :
93 : realT m_v_wind{ 0 }; ///< \f$ C_n^2 \f$ averaged windspeed
94 :
95 : realT m_dir_wind{ 0 }; ///< \f$ C_n^2 \f$ averaged direction
96 :
97 : bool m_z_mean_updated{ false }; ///< whether or not m_z_mean has been updated after changes
98 :
99 : realT m_z_mean{ 0 }; ///< \f$ C_n^2 \f$ averaged layer height
100 :
101 : public:
102 : /** \name PSD Parameters
103 : * @{
104 : */
105 :
106 : /// Get the value of Fried's parameter r_0 at the reference wavelength lam_0.
107 : /**
108 : * \returns the curret value of m_r_0, in m.
109 : */
110 : realT r_0();
111 :
112 : /// Get the value of Fried's parameter r_0 at the specified wavelength.
113 : /**
114 : * \note This does not set the value of r_0.
115 : *
116 : * \returns the current value of m_r_0 * pow(lam, 6/5), in m.
117 : */
118 : realT r_0( const realT &lam /**< [in] the wavelength, in m, at which to calculate r_0*/ );
119 :
120 : /// Set the value of Fried's parameter and the reference wavelength.
121 : /** If the provided reference wavelength is \<=0, then 0.5 microns is used.
122 : *
123 : */
124 : void r_0( const realT &r0, ///< [in] is the new value of r_0, m
125 : const realT &l0 ///< [in] is the new value of lam_0 (in m), if 0 then 0.5e-6 is the default.
126 : );
127 :
128 : /// Get the current value of the reference wavelength.
129 : /** This is the wavelength at which r_0 is specified.
130 : *
131 : * \returns the current value of m_lam_0, in m.
132 : */
133 : realT lam_0();
134 :
135 : /// Get the strength of a single layer.
136 : /**
137 : * \returns the value of m_layer_Cn2[n].
138 : */
139 : realT layer_Cn2( const int n /**< [in] specifies the layer. */ );
140 :
141 : /// Get the vector of layer strengths.
142 : /**
143 : * \returns a copy of the vector of layer strengths: m_layer_Cn2.
144 : */
145 : std::vector<realT> layer_Cn2();
146 :
147 : /// Set the vector layer strengths, possibly calculating r_0.
148 : /**
149 : * If a reference wavelength is specified (l0 > 0), then r_0 is set from the layer strengths
150 : * according to
151 : *
152 : * Regardless of what units the strengths are specified in, they are stored normalized so that
153 : * \f$ \sum_n C_n^2 = 1 \f$.
154 : *
155 : */
156 : error_t layer_Cn2( const std::vector<realT> &cn2, ///< [in] vector containing the layer strengths
157 : const realT l0 = 0 ///< [in] reference wavelength; if positive, also calculate r_0
158 : );
159 :
160 : /// Get the value of the outer scale for a single layer.
161 : /**
162 : * \returns the current value of m_L_0[n], in m.
163 : */
164 : realT L_0( const size_t &n );
165 :
166 : /// Set the vector of layer outer scales.
167 : /**
168 : */
169 : void L_0( const std::vector<realT> &L0 /**< [in] is the new vector of outer scales, in m */ );
170 :
171 : /// Get the vector of outer scales.
172 : /**
173 : * \returns a copy of the vector of layer outer scales, in m
174 : */
175 : std::vector<realT> L_0();
176 :
177 : /// Get the value of the inner scale for a single layer.
178 : /**
179 : * \returns the current value of m_l_0[n], in m.
180 : */
181 : realT l_0( const size_t &n );
182 :
183 : /// Set the vector of layer inner scales.
184 : /**
185 : */
186 : void l_0( const std::vector<realT> &l0 /**< [in] is the new vector of inner scales, in m*/ );
187 :
188 : /// Get the vector of inner scales.
189 : /**
190 : * \returns a copy of the vector of layer inner scales, in m
191 : */
192 : std::vector<realT> l_0();
193 :
194 : /// Set the value of m_nonKolmogorov
195 : /** This flag indicates if non-Kolmogorov turbulence is being modeled.
196 : */
197 : void nonKolmogorov( const bool &nk /**< [in] the value of m_nonKolmogorov*/ );
198 :
199 : /// Return the value of m_nonKolmogorov
200 : /** This flag indicates if non-Kolmogorov turbulence is being modeled.
201 : *
202 : * \returns the current value of m_nonKolmogorov
203 : */
204 : bool nonKolmogorov();
205 :
206 : /// Return the PSD index for a single layer.
207 : /** Satifies the requirements of psdParamsT.
208 : *
209 : * If m_nonKolmogorov is false, this returns
210 : * \f[
211 : \alpha = \frac{11}{3}
212 : \f]
213 : * Otherwise it returns the current value of m_alpha[n].
214 : *
215 : * \returns the PSD index.
216 : */
217 : realT alpha( const size_t &n );
218 :
219 : /// Set the vector of layer PSD indices.
220 : /**
221 : */
222 : void alpha( const std::vector<realT> &alph /**< [in] is the new vector of PSD indices*/ );
223 :
224 : /// Get the vector of PSD indices.
225 : /**
226 : * \returns a copy of the vector of PSD indices
227 : */
228 : std::vector<realT> alpha();
229 :
230 : /// Return the PSD normalization for a single layer.
231 : /** Satifies the requirements of psdParamsT.
232 : *
233 : * If m_nonKolmogorov is false, this returns
234 : * \f[
235 : \beta = \frac{0.0218}{r_0^{5/3}}
236 : \f]
237 : * Otherwise it returns the current value of m_beta[n].
238 : *
239 : * \returns the PSD normalization.
240 : */
241 : realT beta( const size_t &n );
242 :
243 : /// Set the vector of layer PSD normalizations.
244 : /**
245 : */
246 : void beta( const std::vector<realT> &bet /**< [in] is the new vector of PSD normalizations*/ );
247 :
248 : /// Get the vector of PSD normalizations.
249 : /**
250 : * \returns a copy of the vector of PSD normalizations
251 : */
252 : std::vector<realT> beta();
253 :
254 : /// Return the PSD constant for a single layer.
255 : /** Satifies the requirements of psdParamsT.
256 : *
257 : * If m_nonKolmogorov is false, this returns
258 : * \f[
259 : \beta_0 = 0
260 : \f]
261 : * Otherwise it returns the current value of m_beta_0[n].
262 : *
263 : * \returns the PSD constant.
264 : */
265 : realT beta_0( const size_t &n );
266 :
267 : /// Set the vector of layer PSD constants.
268 : /**
269 : */
270 : void beta_0( const std::vector<realT> &bet /**< [in] is the new vector of PSD constants*/ );
271 :
272 : /// Get the vector of PSD constants.
273 : /**
274 : * \returns a copy of the vector of PSD constants
275 : */
276 : std::vector<realT> beta_0();
277 :
278 : /// Get the number of layers
279 : /**
280 : * \returns the size of the m_layer_Cn2 vector if m_nonKolmogrov==false
281 : * \returns the size of the m_alpha vector if m_nonKolmogrov==true
282 : */
283 : size_t n_layers();
284 :
285 : size_t currentLayer();
286 :
287 : void currentLayer( size_t cl );
288 :
289 : /// @}
290 :
291 : /// Get the height of a single layer.
292 : /**
293 : * \returns the height of layer n, in m.
294 : */
295 : realT layer_z( const size_t n /**< [in] specifies the layer. */ );
296 :
297 : /// Get the vector layer heights.
298 : /**
299 : * \returns a copy of the vector of layer heights:m_layer_z, in m.
300 : */
301 : std::vector<realT> layer_z();
302 :
303 : /// Set the vector of layer heights.
304 : /**
305 : */
306 : void layer_z( const std::vector<realT> &layz /**< [in] new vector of layer heights, in m. */ );
307 :
308 : /// Get the height of the observatory.
309 : /**
310 : * \return the current value of m_h_obs, in m.
311 : */
312 : realT h_obs();
313 :
314 : /// Set the height of the observatory.
315 : /**
316 : */
317 : void h_obs( realT nh /**< [in] the new height of the observatory [m] */ );
318 :
319 : /// Get the atmospheric scale height.
320 : /**
321 : * \returns the current value of m_H, in m.
322 : */
323 : realT H();
324 :
325 : /// Set the atmospheric scale height.
326 : /**
327 : */
328 : void H( realT nH /**< [in] the new value of m_H [m] */ );
329 :
330 : /// Get the wind speed of a single layer.
331 : /**
332 : * \returns the value of m_layer_v_wind[n].
333 : */
334 : realT layer_v_wind( const int n /**< [in] specifies the layer. */ );
335 :
336 : /// Get the vector of layer windspeeds
337 : /**
338 : * \returns a copy of the vector of layer windspeeds: m_layer_v_wind.
339 : */
340 : std::vector<realT> layer_v_wind();
341 :
342 : /// Set the vector of layer windspeeds.
343 : void layer_v_wind( const std::vector<realT> &spd /**< [in] the new vector, which is copied to m_layer_v_wind. */ );
344 :
345 : /// Get the wind direction of a single layer.
346 : /**
347 : * \returns the value of m_layer_dir[n], which is the wind direction in that layer in radians.
348 : */
349 : realT layer_dir( const int n /**< [in] specifies the layer. */ );
350 :
351 : /// Get the vector of layer wind directions
352 : /**
353 : * \returns a copy of the vector of layer wind directions: m_layer_dir.
354 : */
355 : std::vector<realT> layer_dir();
356 :
357 : /// Set the vector of layer wind directions.
358 : void layer_dir( const std::vector<realT>
359 : &d /**< [in] the new vector of wind directions in radians, which is copied to m_layer_dir. */ );
360 :
361 : /// Get the 5/3 moment weighted mean wind speed
362 : /** Returns the weighted mean wind speed according to the 5/3's turbulence moment. This is defined as
363 : *
364 : \f[
365 : \bar{v} = \left[\sum_i C_N^2(z_i) v_i^{5/3} \right]^{3/5}
366 : \f]
367 : * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
368 : *
369 : * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
370 : * returns the value of m_v_wind.
371 : *
372 : * \returns the current value of m_v_wind.
373 : */
374 : realT v_wind();
375 :
376 : /// Get the mean wind speed
377 : /** Returns the weighted mean wind speed according to the 5/3's turbulence moment. This is defined as
378 : *
379 : \f[
380 : \bar{v} = \sum_i C_N^2(z_i) v_i
381 : \f]
382 : * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
383 : *
384 : * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
385 : * returns the value of m_v_wind.
386 : *
387 : * \returns the current value of m_v_wind.
388 : */
389 : realT v_wind_mean();
390 :
391 : /// Get the mean-squared wind speed
392 : realT v_wind_mean2();
393 :
394 0 : realT v_max()
395 : {
396 0 : auto max = std::max_element( std::begin( m_layer_v_wind ), std::end( m_layer_v_wind ) );
397 :
398 0 : return *max;
399 : }
400 :
401 : /// Get the weighted mean wind direction
402 : /** Returns the weighted mean wind direction according to the 5/3's turbulence moment. This is defined as
403 : *
404 : \f[
405 : \bar{\theta} = \left[\sum_i C_N^2(z_i) \theta_i^{5/3} \right]^{3/5}
406 : \f]
407 : * See Hardy (1998) Section 3.3.6. \cite hardy_1998.
408 : *
409 : * This is only re-calculated if either m_layer_Cn2 or m_layer_v_wind is changed, otherwise this just
410 : * returns the value of m_dir_wind.
411 : *
412 : * \returns the current value of m_dir_wind.
413 : */
414 : realT dir_wind();
415 :
416 : protected:
417 : /// Recalculate m_v_wind
418 : /** Called by v_wind() whenever m_v_wind_updated is true.
419 : *
420 : * \todo handle dir_wind averaging across 0/2pi
421 : */
422 : void update_v_wind();
423 :
424 : public:
425 : /// Set the weighted mean m_v_wind and renormalize the layer wind speeds
426 : /** Calling this function changes the values of m_layer_v_wind so that the
427 : * Layer averaged 5/3 \f$C_n^2\f$ moment of wind speed is the new value specified by vw.
428 : *
429 : */
430 : void v_wind( const realT &vw /**< [in] the new value of m_v_wind. */ );
431 :
432 : /// Get the weighted mean layer height
433 : /** Returns the weighted layer height according to the 5/3's turbulence moment. This is defined as
434 : \f[
435 : \bar{z} = \left[\sum_i C_N^2(z_i) z_i^{5/3} \right]^{3/5}
436 : \f]
437 : * See Hardy (1998) Section 3.3.6 and 3.7.2. \cite hardy_1998.
438 : *
439 : * This is only re-calculated if either m_layer_Cn2 orm_layer_z is changed, otherwise this just
440 : * returns the value of m_z_mean.
441 : *
442 : * \returns the current value of m_z_mean.
443 : */
444 : realT z_mean();
445 :
446 : protected:
447 : /// Recalculate m_z_mean
448 : /** Called by z_mean() whenever m_z_mean_updated is true.
449 : */
450 : void update_z_mean();
451 :
452 : public:
453 : /// Set the weighted mean m_z_mean and renormalize the layer heights
454 : /** Calling this function changes the values ofm_layer_z so that the
455 : * Layer averaged 5/3 \f$C_n^2\f$ moment of the height is the new value specified.
456 : *
457 : */
458 : void z_mean( const realT &zm /**< [in] is the new value of m_v_wind. */ );
459 :
460 : /// The fraction of the turbulence PSD in phase after Fresnel propagation.
461 : /** See Equation (14) of Guyon (2005) \cite guyon_2005.
462 : *
463 : * \returns the value of the X function.
464 : */
465 : realT X( realT k, ///< [in] the spatial frequency, in inverse meters.
466 : realT lam_sci, ///< [in] is the science observation wavelength.
467 : realT secZ ///< [in] is the secant of the zenith distance.
468 : );
469 :
470 : /// The differential fraction of the turbulence PSD in phase after Fresnel propagation.
471 : /** See Equation (25) of Guyon (2005) \cite guyon_2005.
472 : *
473 : * \returns the value of the dX function.
474 : */
475 : realT dX( realT k, ///< [in] the spatial frequency, in inverse meters.
476 : realT lam_sci, ///< [in] is the science observation wavelength.
477 : realT lam_wfs ///< [in] is the wavefront sensor wavelength.
478 : );
479 :
480 : /// The fraction of the turbulence PSD in amplitude after Fresnel propagation.
481 : /** See Equation (15) of Guyon (2005) \cite guyon_2005.
482 : *
483 : * \returns the value of the Y function.
484 : */
485 : realT Y( realT k, ///< [in] the spatial frequency, in inverse meters.
486 : realT lam_sci, ///< [in] is the science observation wavelength.
487 : realT secZ ///< [in] is the secant of the zenith distance.
488 : );
489 :
490 : /// The differential fraction of the turbulence PSD in amplitude after Fresnel propagation.
491 : /** See Equation (27) of Guyon (2005) \cite guyon_2005.
492 : *
493 : * \returns the value of the dY function.
494 : */
495 : realT dY( realT k, ///< [in] the spatial frequency, in inverse meters.
496 : realT lam_sci, ///< [in] is the science observation wavelength.
497 : realT lam_wfs ///< [in] is the wavefront sensor wavelength.
498 : );
499 :
500 : realT n_air( realT lam /**< [in]The wavelength*/ );
501 :
502 : realT X_Z( realT k, ///< [in] the spatial frequency, in inverse meters
503 : realT lambda_sci, ///< [in] is the science observation wavelength.
504 : realT lambda_wfs, ///< [in] is the wavefront sensor wavelength.
505 : realT secZ ///< [in] is the secant of the zenith distance.
506 : );
507 :
508 : /// Calculate the full-width at half-maximum of a seeing limited image for this atmosphere for a small telescope
509 : /// (ignoring L_0)
510 : /** Calculate the FWHM of a seeing limited image with the current parameters according to Floyd et al. (2010)
511 : \cite floyd_2010 \f[ \epsilon_0 = 0.98\frac{\lambda_{sci}}{r_0(\lambda_sci)}. \f]
512 : *
513 : * \returns the value of the FWHM for the current atmosphere parameters.
514 : */
515 : realT fwhm0( realT lam_sci /**< [in] the wavelength of the science observation. */ );
516 :
517 : /// Calculate the full-width at half-maximum of a seeing limited image for this atmosphere for a large telescope
518 : /// (including L_0)
519 : /** Calculate the FWHM of a seeing limited image with the current parameters according to Floyd et al. (2010)
520 : \cite floyd_2010 \f[ \epsilon_0 = 0.98\frac{\lambda_{sci}}{r_0(\lambda_sci)}. \f]
521 : * If there is an outer scale (_L_0 > 0), then a correction is applied according to Tokovinin (2002)
522 : \cite tokovinin_2002 \f[ \left( \frac{\epsilon_{vK}}{\epsilon_0}\right)^2 = 1 - 2.183\left(
523 : \frac{r_0(\lambda_{sci}}{L_0}\right)^{0.356} \f]
524 : *
525 : *
526 : * \returns the value of the FWHM (\f$ \epsilon_{0/vK} \f$) for the current atmosphere parameters.
527 : */
528 : realT fwhm( realT lam_sci /**< [in] the wavelength of the science observation. */ );
529 :
530 : /// Get the greenwood frequency at the reference wavelength
531 : /**
532 : * \todo derive full value of the constant
533 : */
534 : realT f_g();
535 :
536 : /// Get the greenwood frequency at a specified wavelength
537 : /**
538 : *
539 : */
540 : realT f_g( realT lam_sci /**< [in] the wavelength of the science observation. */ );
541 :
542 : /// Get tau_0 at the reference wavelength
543 : /**
544 : * \todo derive full value of the constant
545 : */
546 : realT tau_0();
547 :
548 : /// Get tau_0 at a specified wavelength.
549 : /**
550 : *
551 : */
552 : realT tau_0( realT lam_sci /**< [in] the wavelength of the science observation. */ );
553 :
554 : /// Scale v_wind so that tau_0 has the specified value at the specified wavelength.
555 : /** Does not modify r_0.
556 : */
557 : void tau_0( realT tau_0, ///< [in] the desired tau_0
558 : realT lam_sci ///< [in] the wavelength of the science observation.
559 : );
560 :
561 : /// Load the default atmosphere model from Guyon (2005).
562 : /** Sets the parameters from Table 4 of Guyon (2005) \cite guyon_2005.
563 : */
564 : void loadGuyon2005();
565 :
566 : /// Load parameters corresponding to the median atmosphere of the GMT site survey at LCO.
567 : /**
568 : */
569 : void loadLCO();
570 :
571 : /// Set a single layer model.
572 : /** Sets all layer vectors to size=1 and populates their fields based on these arguments.
573 : *
574 : */
575 : void setSingleLayer( realT r0, ///< [in] is the new value of r_0
576 : realT lam0, ///< [in] is the new value of lam_0, if 0 then 0.5 microns is the default.
577 : realT L0, ///< [in] the new outer scale
578 : realT l0, ///< [in] the new inner scale
579 : realT lz, ///< [in] the layer height
580 : realT vw, ///< [in] the layer wind-speed
581 : realT dir ///< [in] the layer wind direction.
582 : );
583 :
584 : /// Output current parameters to a stream
585 : /** Prints a formatted list of all current parameters.
586 : *
587 : * \tparam iosT is a std::ostream-like type.
588 : *
589 : * \todo update for new vector components (L_0, etc.)
590 : */
591 : template <typename iosT>
592 : iosT &dumpAtmosphere( iosT &ios /**< [in] a std::ostream-like stream. */ );
593 :
594 : /// \name mx::application support
595 : /** @{
596 : */
597 :
598 : /// Setup the configurator to configure this class.
599 : void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
600 :
601 : /// Load the configuration of this class from a configurator.
602 : /** The complete atmosphere is validated before derived rescalings are applied and again before returning.
603 : *
604 : * \returns `error_t::noerror` for a valid loaded atmosphere, or a typed configuration error.
605 : */
606 : error_t loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
607 :
608 : /// @}
609 : };
610 :
611 : template <typename realT>
612 595 : aoAtmosphere<realT>::aoAtmosphere()
613 : {
614 85 : }
615 :
616 : template <typename realT>
617 79 : error_t aoAtmosphere<realT>::validate() const
618 : {
619 79 : const size_t layerCount = m_layer_Cn2.size();
620 79 : if( layerCount == 0 )
621 : {
622 33 : return internal::mxlib_error_report( error_t::sizeerr, "atmosphere must contain at least one layer" );
623 : }
624 :
625 198 : if( m_L_0.size() != layerCount || m_l_0.size() != layerCount || m_layer_z.size() != layerCount ||
626 198 : m_layer_v_wind.size() != layerCount || m_layer_dir.size() != layerCount )
627 : {
628 21 : return internal::mxlib_error_report( error_t::sizeerr, "atmosphere layer-vector sizes do not match" );
629 : }
630 :
631 90 : if( m_nonKolmogorov &&
632 29 : ( m_alpha.size() != layerCount || m_beta.size() != layerCount || m_beta_0.size() != layerCount ) )
633 : {
634 9 : return internal::mxlib_error_report( error_t::sizeerr, "non-Kolmogorov atmosphere vector sizes do not match" );
635 : }
636 :
637 114 : if( !math::isFinite( m_h_obs ) || m_h_obs < 0 || !math::isFinite( m_H ) || m_H <= 0 ||
638 56 : ( !m_nonKolmogorov &&
639 44 : ( !math::isFinite( m_r_0 ) || m_r_0 <= 0 || !math::isFinite( m_lam_0 ) || m_lam_0 <= 0 ) ) )
640 : {
641 9 : return internal::mxlib_error_report( error_t::invalidconfig, "atmosphere scalar parameters are invalid" );
642 : }
643 :
644 55 : realT totalStrength = 0;
645 55 : bool hasPositiveStrength = false;
646 134 : for( size_t index = 0; index < layerCount; ++index )
647 : {
648 170 : if( !math::isFinite( m_layer_Cn2[index] ) || m_layer_Cn2[index] < 0 || !math::isFinite( m_L_0[index] ) ||
649 84 : !math::isFinite( m_l_0[index] ) || m_l_0[index] < 0 || !math::isFinite( m_layer_z[index] ) ||
650 251 : m_layer_z[index] < 0 || !math::isFinite( m_layer_v_wind[index] ) || m_layer_v_wind[index] < 0 ||
651 81 : !math::isFinite( m_layer_dir[index] ) )
652 : {
653 15 : return internal::mxlib_error_report( error_t::invalidconfig, "atmosphere layer parameters are invalid" );
654 : }
655 :
656 80 : totalStrength += m_layer_Cn2[index];
657 80 : hasPositiveStrength = hasPositiveStrength || m_layer_Cn2[index] > 0;
658 :
659 92 : if( m_nonKolmogorov && ( !math::isFinite( m_alpha[index] ) || !math::isFinite( m_beta[index] ) ||
660 12 : m_beta[index] <= 0 || !math::isFinite( m_beta_0[index] ) || m_beta_0[index] < 0 ) )
661 : {
662 2 : return internal::mxlib_error_report( error_t::invalidconfig,
663 3 : "non-Kolmogorov atmosphere parameters are invalid" );
664 : }
665 : }
666 :
667 49 : if( !hasPositiveStrength || !math::isFinite( totalStrength ) || totalStrength <= 0 )
668 : {
669 0 : return internal::mxlib_error_report( error_t::invalidconfig,
670 0 : "atmosphere must contain a finite positive layer-strength sum" );
671 : }
672 :
673 49 : return error_t::noerror;
674 : }
675 :
676 : template <typename realT>
677 22 : realT aoAtmosphere<realT>::r_0()
678 : {
679 22 : return m_r_0;
680 : }
681 :
682 : template <typename realT>
683 7 : realT aoAtmosphere<realT>::r_0( const realT &lam )
684 : {
685 7 : return m_r_0 * pow( lam / m_lam_0, math::six_fifths<realT>() );
686 : }
687 :
688 : template <typename realT>
689 78 : void aoAtmosphere<realT>::r_0( const realT &r0, const realT &l0 )
690 : {
691 78 : m_r_0 = r0;
692 :
693 78 : if( l0 > 0 )
694 : {
695 78 : m_lam_0 = l0;
696 : }
697 : else
698 : {
699 0 : m_lam_0 = 0.5e-6;
700 : }
701 78 : }
702 :
703 : template <typename realT>
704 4 : realT aoAtmosphere<realT>::lam_0()
705 : {
706 4 : return m_lam_0;
707 : }
708 :
709 : template <typename realT>
710 3 : realT aoAtmosphere<realT>::layer_Cn2( const int n )
711 : {
712 3 : return m_layer_Cn2[n];
713 : }
714 :
715 : template <typename realT>
716 14 : std::vector<realT> aoAtmosphere<realT>::layer_Cn2()
717 : {
718 14 : return m_layer_Cn2;
719 : }
720 :
721 : template <typename realT>
722 84 : error_t aoAtmosphere<realT>::layer_Cn2( const std::vector<realT> &cn2, const realT l0 )
723 : {
724 84 : if( cn2.empty() || !math::isFinite( l0 ) || l0 < 0 )
725 : {
726 4 : return internal::mxlib_error_report( error_t::invalidarg,
727 6 : "layer strengths must be nonempty and reference wavelength nonnegative" );
728 : }
729 :
730 82 : realT layer_norm = 0;
731 183 : for( size_t i = 0; i < cn2.size(); ++i )
732 : {
733 103 : if( !math::isFinite( cn2[i] ) || cn2[i] < 0 )
734 : {
735 4 : return internal::mxlib_error_report( error_t::invalidarg,
736 6 : "layer strengths must be finite and nonnegative" );
737 : }
738 101 : layer_norm += cn2[i];
739 : }
740 :
741 80 : if( !math::isFinite( layer_norm ) || layer_norm <= 0 )
742 : {
743 9 : return internal::mxlib_error_report( error_t::invalidarg, "layer strengths must have a finite positive sum" );
744 : }
745 :
746 77 : std::vector<realT> normalized = cn2;
747 170 : for( size_t i = 0; i < normalized.size(); ++i )
748 : {
749 93 : normalized[i] /= layer_norm;
750 : }
751 :
752 77 : if( l0 > 0 )
753 : {
754 0 : const realT r0 = 1.0 / pow( layer_norm * 5.520e13, math::three_fifths<realT>() );
755 0 : if( !math::isFinite( r0 ) || r0 <= 0 )
756 : {
757 0 : return internal::mxlib_error_report( error_t::invalidarg,
758 0 : "layer strengths produce an invalid Fried parameter" );
759 : }
760 0 : m_r_0 = r0;
761 0 : m_lam_0 = l0;
762 : }
763 :
764 77 : m_layer_Cn2 = std::move( normalized );
765 :
766 77 : m_v_wind_updated = false;
767 77 : m_z_mean_updated = false;
768 :
769 77 : return error_t::noerror;
770 77 : }
771 :
772 : template <typename realT>
773 288822 : realT aoAtmosphere<realT>::L_0( const size_t &n )
774 : {
775 288822 : return m_L_0[n];
776 : }
777 :
778 : template <typename realT>
779 78 : void aoAtmosphere<realT>::L_0( const std::vector<realT> &L_0 )
780 : {
781 78 : m_L_0 = L_0;
782 78 : }
783 :
784 : template <typename realT>
785 1 : std::vector<realT> aoAtmosphere<realT>::L_0()
786 : {
787 1 : return m_L_0;
788 : }
789 :
790 : template <typename realT>
791 5 : realT aoAtmosphere<realT>::l_0( const size_t &n )
792 : {
793 5 : return m_l_0[n];
794 : }
795 :
796 : template <typename realT>
797 77 : void aoAtmosphere<realT>::l_0( const std::vector<realT> &l_0 )
798 : {
799 77 : m_l_0 = l_0;
800 77 : }
801 :
802 : template <typename realT>
803 1 : std::vector<realT> aoAtmosphere<realT>::l_0()
804 : {
805 1 : return m_l_0;
806 : }
807 :
808 : template <typename realT>
809 6 : void aoAtmosphere<realT>::nonKolmogorov( const bool &nk )
810 : {
811 6 : m_nonKolmogorov = nk;
812 6 : }
813 :
814 : template <typename realT>
815 15 : bool aoAtmosphere<realT>::nonKolmogorov()
816 : {
817 15 : return m_nonKolmogorov;
818 : }
819 :
820 : template <typename realT>
821 96616 : realT aoAtmosphere<realT>::alpha( const size_t &n )
822 : {
823 96616 : if( !m_nonKolmogorov )
824 : {
825 96615 : return math::eleven_thirds<realT>();
826 : }
827 : else
828 : {
829 1 : return m_alpha[n];
830 : }
831 : }
832 :
833 : template <typename realT>
834 2 : void aoAtmosphere<realT>::alpha( const std::vector<realT> &alph )
835 : {
836 2 : m_nonKolmogorov = true;
837 2 : m_alpha = alph;
838 2 : }
839 :
840 : template <typename realT>
841 0 : std::vector<realT> aoAtmosphere<realT>::alpha()
842 : {
843 0 : return m_alpha;
844 : }
845 :
846 : template <typename realT>
847 96451 : realT aoAtmosphere<realT>::beta( const size_t &n )
848 : {
849 96451 : if( !m_nonKolmogorov )
850 : {
851 96450 : return constants::a_PSD<realT>() * pow( m_r_0, -math::five_thirds<realT>() );
852 : ;
853 : }
854 : else
855 : {
856 1 : return m_beta[n];
857 : }
858 : }
859 :
860 : template <typename realT>
861 3 : void aoAtmosphere<realT>::beta( const std::vector<realT> &bet )
862 : {
863 3 : m_nonKolmogorov = true;
864 3 : m_beta = bet;
865 3 : }
866 :
867 : template <typename realT>
868 0 : std::vector<realT> aoAtmosphere<realT>::beta()
869 : {
870 0 : return m_beta;
871 : }
872 :
873 : template <typename realT>
874 96451 : realT aoAtmosphere<realT>::beta_0( const size_t &n )
875 : {
876 96451 : if( !m_nonKolmogorov )
877 : {
878 96450 : return 0;
879 : }
880 : else
881 : {
882 1 : return m_beta_0[n];
883 : }
884 : }
885 :
886 : template <typename realT>
887 2 : void aoAtmosphere<realT>::beta_0( const std::vector<realT> &bet )
888 : {
889 2 : m_nonKolmogorov = true;
890 2 : m_beta_0 = bet;
891 2 : }
892 :
893 : template <typename realT>
894 0 : std::vector<realT> aoAtmosphere<realT>::beta_0()
895 : {
896 0 : return m_beta_0;
897 : }
898 :
899 : template <typename realT>
900 49 : size_t aoAtmosphere<realT>::n_layers()
901 : {
902 49 : if( m_nonKolmogorov )
903 0 : return m_alpha.size();
904 : else
905 49 : return m_layer_Cn2.size();
906 : }
907 :
908 : template <typename realT>
909 3 : realT aoAtmosphere<realT>::layer_z( const size_t n )
910 : {
911 3 : return m_layer_z[n];
912 : }
913 :
914 : template <typename realT>
915 78 : void aoAtmosphere<realT>::layer_z( const std::vector<realT> &z )
916 : {
917 78 : m_layer_z = z;
918 78 : m_z_mean_updated = false;
919 78 : }
920 :
921 : template <typename realT>
922 1 : std::vector<realT> aoAtmosphere<realT>::layer_z()
923 : {
924 1 : return m_layer_z;
925 : }
926 :
927 : template <typename realT>
928 2 : realT aoAtmosphere<realT>::h_obs()
929 : {
930 2 : return m_h_obs;
931 : }
932 :
933 : template <typename realT>
934 3 : void aoAtmosphere<realT>::h_obs( realT nh )
935 : {
936 3 : m_h_obs = nh;
937 3 : }
938 :
939 : template <typename realT>
940 2 : realT aoAtmosphere<realT>::H()
941 : {
942 2 : return m_H;
943 : }
944 :
945 : template <typename realT>
946 1 : void aoAtmosphere<realT>::H( realT nH )
947 : {
948 1 : m_H = nH;
949 1 : }
950 :
951 : template <typename realT>
952 96501 : realT aoAtmosphere<realT>::layer_v_wind( const int n )
953 : {
954 96501 : return m_layer_v_wind[n];
955 : }
956 :
957 : template <typename realT>
958 1 : std::vector<realT> aoAtmosphere<realT>::layer_v_wind()
959 : {
960 1 : return m_layer_v_wind;
961 : }
962 :
963 : template <typename realT>
964 80 : void aoAtmosphere<realT>::layer_v_wind( const std::vector<realT> &v )
965 : {
966 80 : m_layer_v_wind = v;
967 80 : m_v_wind_updated = false;
968 80 : }
969 :
970 : template <typename realT>
971 24 : realT aoAtmosphere<realT>::layer_dir( const int n )
972 : {
973 24 : return m_layer_dir[n];
974 : }
975 :
976 : template <typename realT>
977 1 : std::vector<realT> aoAtmosphere<realT>::layer_dir()
978 : {
979 1 : return m_layer_dir;
980 : }
981 :
982 : template <typename realT>
983 79 : void aoAtmosphere<realT>::layer_dir( const std::vector<realT> &dir )
984 : {
985 79 : m_layer_dir = dir;
986 79 : m_v_wind_updated = false;
987 79 : }
988 :
989 : template <typename realT>
990 5 : realT aoAtmosphere<realT>::v_wind()
991 : {
992 5 : if( m_v_wind_updated == false )
993 3 : update_v_wind();
994 5 : return m_v_wind;
995 : }
996 :
997 : template <typename realT>
998 1 : realT aoAtmosphere<realT>::dir_wind()
999 : {
1000 1 : if( m_v_wind_updated == false )
1001 0 : update_v_wind();
1002 1 : return m_dir_wind;
1003 : }
1004 :
1005 : template <typename realT>
1006 3 : void aoAtmosphere<realT>::update_v_wind()
1007 : {
1008 :
1009 3 : if( m_layer_v_wind.size() == 0 || m_layer_Cn2.size() == 0 )
1010 : {
1011 0 : m_v_wind = 0;
1012 0 : m_dir_wind = 0;
1013 0 : m_v_wind_updated = true;
1014 0 : return;
1015 : }
1016 :
1017 3 : m_v_wind = 0;
1018 :
1019 3 : realT s = 0;
1020 3 : realT c = 0;
1021 :
1022 10 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1023 : {
1024 7 : m_v_wind += m_layer_Cn2[i] * pow( m_layer_v_wind[i], math::five_thirds<realT>() );
1025 7 : s += pow( m_layer_v_wind[i], math::five_thirds<realT>() ) *
1026 7 : sin( m_layer_dir[i] ); // pow( sin(m_layer_dir[i]), math::five_thirds<realT>() );
1027 7 : c += pow( m_layer_v_wind[i], math::five_thirds<realT>() ) *
1028 7 : cos( m_layer_dir[i] ); // pow( cos(m_layer_dir[i]), math::five_thirds<realT>() );
1029 : }
1030 :
1031 3 : m_v_wind = pow( m_v_wind, math::three_fifths<realT>() );
1032 :
1033 : // m_dir_wind = atan(pow(s, math::three_fifths<realT>()) / pow(c, math::three_fifths<realT>()));
1034 3 : m_dir_wind = atan( s / c );
1035 3 : if( m_dir_wind < 0 )
1036 0 : m_dir_wind += math::pi<realT>();
1037 :
1038 3 : m_v_wind_updated = true;
1039 : }
1040 :
1041 : template <typename realT>
1042 1 : void aoAtmosphere<realT>::v_wind( const realT &vw )
1043 : {
1044 1 : if( m_v_wind_updated == false )
1045 0 : update_v_wind();
1046 :
1047 1 : realT vw_old = m_v_wind;
1048 :
1049 1 : m_v_wind = vw;
1050 :
1051 : // Now update the layers if needed
1052 1 : if( m_layer_v_wind.size() > 0 )
1053 : {
1054 4 : for( size_t i = 0; i < m_layer_v_wind.size(); ++i )
1055 : {
1056 3 : m_layer_v_wind[i] = m_layer_v_wind[i] * ( vw / vw_old );
1057 : }
1058 : }
1059 1 : }
1060 :
1061 : template <typename realT>
1062 4 : realT aoAtmosphere<realT>::z_mean()
1063 : {
1064 4 : if( m_z_mean_updated == false )
1065 3 : update_z_mean();
1066 4 : return m_z_mean;
1067 : }
1068 :
1069 : template <typename realT>
1070 3 : void aoAtmosphere<realT>::update_z_mean()
1071 : {
1072 3 : if( m_layer_z.size() == 0 || m_layer_Cn2.size() == 0 )
1073 : {
1074 0 : m_z_mean = 0;
1075 0 : m_z_mean_updated = true;
1076 0 : return;
1077 : }
1078 :
1079 3 : m_z_mean = 0;
1080 :
1081 10 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1082 : {
1083 7 : m_z_mean += m_layer_Cn2[i] * pow( m_layer_z[i], math::five_thirds<realT>() );
1084 : }
1085 :
1086 3 : m_z_mean = pow( m_z_mean, math::three_fifths<realT>() );
1087 :
1088 3 : m_z_mean_updated = true;
1089 : }
1090 :
1091 : template <typename realT>
1092 1 : void aoAtmosphere<realT>::z_mean( const realT &zm )
1093 : {
1094 1 : if( m_z_mean_updated == false )
1095 0 : update_z_mean();
1096 :
1097 1 : realT zh_old = m_z_mean;
1098 :
1099 1 : m_z_mean = zm;
1100 :
1101 : // Now update the layers if needed
1102 1 : if( m_layer_z.size() > 0 )
1103 : {
1104 4 : for( size_t i = 0; i < m_layer_z.size(); ++i )
1105 : {
1106 3 : m_layer_z[i] = m_layer_z[i] * ( zm / zh_old );
1107 : }
1108 : }
1109 1 : }
1110 :
1111 : template <typename realT>
1112 0 : realT aoAtmosphere<realT>::X( realT k, realT lam_sci, realT secZ )
1113 : {
1114 0 : realT c = 0;
1115 :
1116 0 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1117 : {
1118 0 : c += m_layer_Cn2[i] * pow( cos( math::pi<realT>() * k * k * lam_sci * m_layer_z[i] * secZ ), 2 );
1119 : }
1120 :
1121 0 : return c;
1122 : }
1123 :
1124 : template <typename realT>
1125 0 : realT aoAtmosphere<realT>::dX( realT f, realT lam_sci, realT lam_wfs )
1126 : {
1127 0 : realT c = 0;
1128 :
1129 0 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1130 : {
1131 0 : c += m_layer_Cn2[i] * pow( ( cos( math::pi<realT>() * f * f * lam_sci * m_layer_z[i] ) -
1132 0 : cos( math::pi<realT>() * f * f * lam_wfs * m_layer_z[i] ) ),
1133 : 2 );
1134 : }
1135 :
1136 0 : return c;
1137 : }
1138 :
1139 : template <typename realT>
1140 0 : realT aoAtmosphere<realT>::Y( realT k, realT lam_sci, realT secZ )
1141 : {
1142 0 : realT c = 0;
1143 :
1144 0 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1145 : {
1146 0 : c += m_layer_Cn2[i] * pow( sin( math::pi<realT>() * k * k * lam_sci * m_layer_z[i] * secZ ), 2 );
1147 : }
1148 0 : return c;
1149 : }
1150 :
1151 : template <typename realT>
1152 0 : realT aoAtmosphere<realT>::dY( realT f, realT lam_sci, realT lam_wfs )
1153 : {
1154 0 : realT c = 0;
1155 :
1156 0 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1157 : {
1158 0 : c += m_layer_Cn2[i] * pow( ( sin( math::pi<realT>() * f * f * lam_sci * m_layer_z[i] ) -
1159 0 : sin( math::pi<realT>() * f * f * lam_wfs * m_layer_z[i] ) ),
1160 : 2 );
1161 : }
1162 :
1163 0 : return c;
1164 : }
1165 :
1166 : template <typename realT>
1167 0 : realT aoAtmosphere<realT>::n_air( realT lambda )
1168 : {
1169 0 : realT ll2 = static_cast<realT>( 1 ) / pow( lambda / 1e-6, 2 );
1170 :
1171 0 : return 1.0 + 8.34213e-5 + 0.0240603 / ( 130.0 - ll2 ) + 0.00015997 / ( 38.9 - ll2 );
1172 : }
1173 :
1174 : template <typename realT>
1175 0 : realT aoAtmosphere<realT>::X_Z( realT k, realT lambda_i, realT lambda_wfs, realT secZ )
1176 : {
1177 0 : realT c = 0;
1178 0 : realT sinZ = sqrt( 1.0 - pow( 1.0 / secZ, 2 ) );
1179 0 : realT tanZ = sinZ * secZ;
1180 0 : realT x0 = ( n_air( lambda_wfs ) - n_air( lambda_i ) ) * m_H * tanZ * secZ;
1181 : realT x;
1182 :
1183 0 : for( size_t i = 0; i < m_layer_Cn2.size(); ++i )
1184 : {
1185 0 : x = x0 * ( 1 - exp( ( m_layer_z[i] + m_h_obs ) / m_H ) );
1186 0 : c += m_layer_Cn2[i] * pow( cos( math::pi<realT>() * k * k * lambda_i * m_layer_z[i] * secZ ), 2 ) *
1187 0 : pow( sin( math::pi<realT>() * x * k * cos( 0. * 3.14 / 180. ) ), 2 );
1188 : }
1189 :
1190 0 : return 4 * c;
1191 : }
1192 :
1193 : template <typename realT>
1194 0 : realT aoAtmosphere<realT>::fwhm0( realT lam_sci )
1195 : {
1196 0 : realT r0lam = r_0( lam_sci );
1197 :
1198 0 : realT fwhm = 0.98 * ( lam_sci / r0lam );
1199 :
1200 0 : return fwhm;
1201 : }
1202 :
1203 : template <typename realT>
1204 1 : realT aoAtmosphere<realT>::fwhm( realT lam_sci )
1205 : {
1206 1 : realT r0lam = r_0( lam_sci );
1207 :
1208 1 : realT fwhm = 0.98 * ( lam_sci / r0lam );
1209 :
1210 : ///\todo this needs to handle layers with different L_0
1211 1 : if( L_0( 0 ) > 0 )
1212 1 : fwhm *= sqrt( 1 - 2.183 * pow( r0lam / L_0( 0 ), 0.356 ) );
1213 :
1214 1 : return fwhm;
1215 : }
1216 :
1217 : template <typename realT>
1218 0 : realT aoAtmosphere<realT>::f_g()
1219 : {
1220 0 : return 0.428 * v_wind() / m_r_0;
1221 : }
1222 :
1223 : template <typename realT>
1224 1 : realT aoAtmosphere<realT>::f_g( realT lam_sci )
1225 : {
1226 1 : return 0.428 * pow( m_lam_0 / lam_sci, math::six_fifths<realT>() ) * v_wind() / m_r_0;
1227 : }
1228 :
1229 : template <typename realT>
1230 0 : realT aoAtmosphere<realT>::tau_0()
1231 : {
1232 0 : return 0.134 / f_g();
1233 : }
1234 :
1235 : template <typename realT>
1236 1 : realT aoAtmosphere<realT>::tau_0( realT lam_sci )
1237 : {
1238 1 : return 0.134 / f_g( lam_sci );
1239 : }
1240 :
1241 : template <typename realT>
1242 0 : void aoAtmosphere<realT>::tau_0( realT tau_0, realT lam_sci )
1243 : {
1244 0 : realT vw = ( 0.134 / tau_0 ) / 0.428 * m_r_0 * pow( lam_sci / m_lam_0, math::six_fifths<realT>() );
1245 0 : v_wind( vw );
1246 0 : }
1247 :
1248 : template <typename realT>
1249 1 : void aoAtmosphere<realT>::loadGuyon2005()
1250 : {
1251 2 : layer_Cn2( { 0.2283, 0.0883, 0.0666, 0.1458, 0.3350, 0.1350 } );
1252 2 : layer_z( { 500, 1000, 2000, 4000, 8000, 16000 } );
1253 2 : layer_v_wind( { 10., 10., 10., 10., 10., 10. } );
1254 2 : layer_dir( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } );
1255 2 : L_0( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } );
1256 2 : l_0( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } );
1257 :
1258 1 : r_0( 0.2, 0.5e-6 );
1259 :
1260 1 : h_obs( 4200 );
1261 1 : }
1262 :
1263 : template <typename realT>
1264 1 : void aoAtmosphere<realT>::loadLCO()
1265 : {
1266 2 : layer_Cn2( { 0.42, 0.029, 0.062, 0.16, 0.11, 0.10, 0.12 } );
1267 2 : layer_z( { 250., 500., 1000., 2000., 4000., 8000., 16000. } );
1268 2 : layer_v_wind( { 10.0, 10.0, 20.0, 20.0, 25.0, 30.0, 25.0 } );
1269 2 : layer_dir( { 1.05, 1.05, 1.31, 1.31, 1.75, 1.92, 1.75 } );
1270 :
1271 1 : r_0( 0.16, 0.5e-6 );
1272 :
1273 2 : L_0( { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 } );
1274 2 : l_0( { 0.0, 0, 0, 0, 0, 0, 0 } );
1275 :
1276 1 : h_obs( 2400 );
1277 1 : }
1278 :
1279 : template <typename realT>
1280 72 : void aoAtmosphere<realT>::setSingleLayer( realT r0, realT lam0, realT L0, realT l0, realT lz, realT vw, realT dir )
1281 : {
1282 72 : r_0( r0, lam0 );
1283 72 : L_0( std::vector<realT>( { L0 } ) );
1284 72 : l_0( std::vector<realT>( { l0 } ) );
1285 72 : layer_Cn2( std::vector<realT>( { 1 } ) );
1286 72 : layer_z( std::vector<realT>( { lz } ) );
1287 72 : layer_v_wind( std::vector<realT>( { vw } ) );
1288 72 : layer_dir( std::vector<realT>( { dir } ) );
1289 72 : }
1290 :
1291 : template <typename realT>
1292 : template <typename iosT>
1293 1 : iosT &aoAtmosphere<realT>::dumpAtmosphere( iosT &ios )
1294 : {
1295 1 : ios << "# Atmosphere Parameters:\n";
1296 1 : ios << "# nonKolmogorov = " << std::boolalpha << nonKolmogorov() << '\n';
1297 1 : ios << "# n_layers = " << n_layers() << '\n';
1298 :
1299 1 : if( !m_nonKolmogorov )
1300 : {
1301 1 : ios << "# r_0 = " << r_0() << '\n';
1302 1 : ios << "# lam_0 = " << lam_0() << '\n';
1303 1 : ios << "# tau_0 = " << tau_0( lam_0() ) << '\n';
1304 1 : ios << "# FWHM = " << fwhm( lam_0() ) << '\n';
1305 1 : ios << "# layer_Cn2 = ";
1306 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1307 0 : ios << layer_Cn2()[i] << ", ";
1308 1 : ios << layer_Cn2()[n_layers() - 1] << '\n';
1309 : }
1310 1 : if( m_nonKolmogorov )
1311 : {
1312 0 : ios << "# alpha = ";
1313 0 : for( size_t i = 0; i < n_layers() - 1; ++i )
1314 0 : ios << alpha()[i] << ", ";
1315 0 : ios << alpha()[n_layers() - 1] << '\n';
1316 0 : ios << "# beta = ";
1317 0 : for( size_t i = 0; i < n_layers() - 1; ++i )
1318 0 : ios << beta()[i] << ", ";
1319 0 : ios << beta()[n_layers() - 1] << '\n';
1320 0 : ios << "# beta_0 = ";
1321 0 : for( size_t i = 0; i < n_layers() - 1; ++i )
1322 0 : ios << beta_0()[i] << ", ";
1323 0 : ios << beta_0()[n_layers() - 1] << '\n';
1324 : }
1325 :
1326 1 : ios << "# L_0 = ";
1327 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1328 0 : ios << L_0()[i] << ", ";
1329 1 : ios << L_0()[n_layers() - 1] << '\n';
1330 1 : ios << "# l_0 = ";
1331 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1332 0 : ios << l_0()[i] << ", ";
1333 1 : ios << l_0()[n_layers() - 1] << '\n';
1334 :
1335 1 : ios << "# layer_v_wind = ";
1336 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1337 0 : ios << layer_v_wind()[i] << ", ";
1338 1 : ios << layer_v_wind()[n_layers() - 1] << '\n';
1339 1 : ios << "# layer_dir = ";
1340 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1341 0 : ios << layer_dir()[i] << ", ";
1342 1 : ios << layer_dir()[n_layers() - 1] << '\n';
1343 1 : ios << "# mean v_wind = " << v_wind() << '\n';
1344 1 : ios << "# mean dir_wind = " << dir_wind() << '\n';
1345 :
1346 1 : if( !m_nonKolmogorov )
1347 : {
1348 1 : ios << "# layer_z = ";
1349 1 : for( size_t i = 0; i < n_layers() - 1; ++i )
1350 0 : ios << layer_z()[i] << ", ";
1351 1 : ios << layer_z()[n_layers() - 1] << '\n';
1352 1 : ios << "# mean z = " << z_mean() << '\n';
1353 1 : ios << "# h_obs = " << h_obs() << '\n';
1354 1 : ios << "# H = " << H() << '\n';
1355 : }
1356 :
1357 1 : return ios;
1358 : }
1359 :
1360 : template <typename realT>
1361 20 : void aoAtmosphere<realT>::setupConfig( app::appConfigurator &config )
1362 : {
1363 : using namespace mx::app;
1364 :
1365 280 : config.add( "atm.r_0", "", "atm.r_0", argType::Required, "atm", "r_0", false, "real", "Fried's parameter [m]" );
1366 20 : config.add( "atm.lam_0",
1367 40 : "",
1368 40 : "atm.lam_0",
1369 : argType::Required,
1370 40 : "atm",
1371 40 : "lam_0",
1372 : false,
1373 40 : "real",
1374 60 : "The reference wavlength for r_0 [m]" );
1375 20 : config.add( "atm.L_0",
1376 40 : "",
1377 40 : "atm.L_0",
1378 : argType::Required,
1379 40 : "atm",
1380 40 : "L_0",
1381 : false,
1382 40 : "vector<real>",
1383 60 : "Layer outer scales [m]" );
1384 20 : config.add( "atm.l_0",
1385 40 : "",
1386 40 : "atm.l_0",
1387 : argType::Required,
1388 40 : "atm",
1389 40 : "l_0",
1390 : false,
1391 40 : "vector<real>",
1392 60 : "Layer inner scales [m]" );
1393 20 : config.add( "atm.layer_z",
1394 40 : "",
1395 40 : "atm.layer_z",
1396 : argType::Required,
1397 40 : "atm",
1398 40 : "layer_z",
1399 : false,
1400 40 : "vector<real>",
1401 60 : "layer heights [m]" );
1402 20 : config.add( "atm.h_obs",
1403 40 : "",
1404 40 : "atm.h_obs",
1405 : argType::Required,
1406 40 : "atm",
1407 40 : "h_obs",
1408 : false,
1409 40 : "real",
1410 60 : "height of observatory [m]" );
1411 280 : config.add( "atm.H", "", "atm.H", argType::Required, "atm", "H", false, "real", "atmospheric scale heights [m]" );
1412 20 : config.add( "atm.layer_Cn2",
1413 40 : "",
1414 40 : "atm.layer_Cn2",
1415 : argType::Required,
1416 40 : "atm",
1417 40 : "layer_Cn2",
1418 : false,
1419 40 : "vector<real>",
1420 60 : "Layer Cn^2. Note that this does not set r_0." );
1421 20 : config.add( "atm.layer_v_wind",
1422 40 : "",
1423 40 : "atm.layer_v_wind",
1424 : argType::Required,
1425 40 : "atm",
1426 40 : "layer_v_wind",
1427 : false,
1428 40 : "vector<real>",
1429 60 : "Layer wind speeds [m/s]" );
1430 20 : config.add( "atm.layer_dir",
1431 40 : "",
1432 40 : "atm.layer_dir",
1433 : argType::Required,
1434 40 : "atm",
1435 40 : "layer_dir",
1436 : false,
1437 40 : "vector<real>",
1438 60 : "Layer wind directions [rad]" );
1439 20 : config.add( "atm.v_wind",
1440 40 : "",
1441 40 : "atm.v_wind",
1442 : argType::Required,
1443 40 : "atm",
1444 40 : "v_wind",
1445 : false,
1446 40 : "real",
1447 60 : "Mean windspeed (5/3 momement), rescales layer speeds [m/s]" );
1448 20 : config.add( "atm.tau_0",
1449 40 : "",
1450 40 : "atm.tau_0",
1451 : argType::Required,
1452 40 : "atm",
1453 40 : "tau_0",
1454 : false,
1455 40 : "real",
1456 60 : "AO time constant, sets v_wind and rescales layer speeds. [s]" );
1457 20 : config.add( "atm.z_mean",
1458 40 : "",
1459 40 : "atm.z_mean",
1460 : argType::Required,
1461 40 : "atm",
1462 40 : "z_mean",
1463 : false,
1464 40 : "real",
1465 60 : "Mean layer height (5/3 momemnt), rescales layer heights [m/s]" );
1466 20 : config.add( "atm.nonKolmogorov",
1467 40 : "",
1468 40 : "atm.nonKolmogorov",
1469 : argType::Required,
1470 40 : "atm",
1471 40 : "nonKolmogorov",
1472 : false,
1473 40 : "bool",
1474 60 : "Set to use a non-Kolmogorov PSD. See alpha and beta." );
1475 20 : config.add( "atm.alpha",
1476 40 : "",
1477 40 : "atm.alpha",
1478 : argType::Required,
1479 40 : "atm",
1480 40 : "alpha",
1481 : false,
1482 40 : "vector<real>",
1483 60 : "Non-kolmogorov PSD exponent." );
1484 20 : config.add( "atm.beta",
1485 40 : "",
1486 40 : "atm.beta",
1487 : argType::Required,
1488 40 : "atm",
1489 40 : "beta",
1490 : false,
1491 40 : "vector<real>",
1492 60 : "Non-kolmogorov PSD normalization." );
1493 20 : config.add( "atm.beta_0",
1494 40 : "",
1495 40 : "atm.beta_0",
1496 : argType::Required,
1497 40 : "atm",
1498 40 : "beta_0",
1499 : false,
1500 40 : "vector<real>",
1501 40 : "Non-kolmogorov PSD constant." );
1502 20 : }
1503 :
1504 : template <typename realT>
1505 20 : error_t aoAtmosphere<realT>::loadConfig( app::appConfigurator &config )
1506 : {
1507 : // Here "has side effecs" means that the set function does more than simply copy the value.
1508 :
1509 : // The order of lam_0, Cn2, and r_0 is so that r_0 overrides the value set with Cn2 if lam_0 != 0.
1510 : // lam_0 comes first because it calibrates r0 and Cn2
1511 :
1512 : // lam_0
1513 20 : config( m_lam_0, "atm.lam_0" );
1514 :
1515 : // layer_Cn2
1516 20 : std::vector<realT> lcn2 = m_layer_Cn2;
1517 40 : config( lcn2, "atm.layer_Cn2" );
1518 40 : if( config.isSet( "atm.layer_Cn2" ) )
1519 : {
1520 3 : const error_t strengthStatus = layer_Cn2( lcn2 );
1521 3 : if( strengthStatus != error_t::noerror )
1522 : {
1523 1 : return strengthStatus;
1524 : }
1525 : }
1526 :
1527 19 : realT r0 = r_0();
1528 38 : config( r0, "atm.r_0" );
1529 38 : if( config.isSet( "atm.r_0" ) )
1530 2 : r_0( r0, m_lam_0 );
1531 :
1532 38 : config( m_L_0, "atm.L_0" );
1533 :
1534 19 : config( m_l_0, "atm.l_0" );
1535 :
1536 : // Has side effects:
1537 19 : std::vector<realT> layz = m_layer_z;
1538 38 : config( layz, "atm.layer_z" ); // Do this no matter what to record source
1539 38 : if( config.isSet( "atm.layer_z" ) )
1540 2 : layer_z( layz ); // but only call this if changed
1541 :
1542 38 : config( m_h_obs, "atm.h_obs" );
1543 19 : config( m_H, "atm.H" );
1544 :
1545 : // Has side effects:
1546 19 : std::vector<realT> lvw = m_layer_v_wind;
1547 38 : config( lvw, "atm.layer_v_wind" ); // Do this no matter what to record source
1548 38 : if( config.isSet( "atm.layer_v_wind" ) )
1549 2 : layer_v_wind( lvw ); // but only call this if changed
1550 :
1551 : // Has side effects:
1552 19 : std::vector<realT> ld = m_layer_dir;
1553 38 : config( ld, "atm.layer_dir" ); // Do this no matter what to record source
1554 38 : if( config.isSet( "atm.layer_dir" ) )
1555 2 : layer_dir( ld ); // but only call this if changed
1556 :
1557 19 : realT vw = 0;
1558 19 : config( vw, "atm.v_wind" ); // Do this no matter what to record source
1559 :
1560 19 : realT t0 = 0;
1561 19 : config( t0, "atm.tau_0" ); // Do this no matter what to record source
1562 19 : realT zm = 0;
1563 38 : config( zm, "atm.z_mean" ); // Do this no matter what to record source
1564 :
1565 19 : config( m_nonKolmogorov, "atm.nonKolmogorov" );
1566 :
1567 19 : std::vector<realT> a = m_alpha;
1568 38 : config( a, "atm.alpha" );
1569 38 : if( config.isSet( "atm.alpha" ) )
1570 1 : alpha( a ); // this sets m_nonKolmogorov
1571 :
1572 19 : std::vector<realT> b = m_beta;
1573 38 : config( b, "atm.beta" );
1574 38 : if( config.isSet( "atm.beta" ) )
1575 1 : beta( b ); // this sets m_nonKolmogorov
1576 :
1577 19 : std::vector<realT> b0 = m_beta_0;
1578 38 : config( b0, "atm.beta_0" );
1579 38 : if( config.isSet( "atm.beta_0" ) )
1580 1 : beta_0( b0 ); // this sets m_nonKolmogorov
1581 :
1582 19 : error_t status = validate();
1583 19 : if( status != error_t::noerror )
1584 : {
1585 9 : return status;
1586 : }
1587 :
1588 20 : if( config.isSet( "atm.v_wind" ) )
1589 : {
1590 1 : if( !math::isFinite( vw ) || vw <= 0 )
1591 : {
1592 0 : return internal::mxlib_error_report( error_t::invalidconfig,
1593 0 : "configured mean wind speed must be finite and positive" );
1594 : }
1595 :
1596 1 : if( v_wind() <= 0 )
1597 : {
1598 0 : return internal::mxlib_error_report( error_t::invalidconfig,
1599 0 : "a static atmosphere cannot be rescaled to positive mean wind" );
1600 : }
1601 1 : v_wind( vw );
1602 : }
1603 :
1604 20 : if( config.isSet( "atm.tau_0" ) )
1605 : {
1606 0 : if( !math::isFinite( t0 ) || t0 <= 0 || !math::isFinite( m_r_0 ) || m_r_0 <= 0 || !math::isFinite( m_lam_0 ) ||
1607 0 : m_lam_0 <= 0 || v_wind() <= 0 )
1608 : {
1609 0 : return internal::mxlib_error_report(
1610 0 : error_t::invalidconfig,
1611 0 : "configured atmosphere time constant requires valid Fried, wavelength, and wind values" );
1612 : }
1613 0 : tau_0( t0, m_lam_0 );
1614 : }
1615 :
1616 20 : if( config.isSet( "atm.z_mean" ) )
1617 : {
1618 1 : if( !math::isFinite( zm ) || zm < 0 )
1619 : {
1620 0 : return internal::mxlib_error_report( error_t::invalidconfig,
1621 0 : "configured mean layer height must be finite and nonnegative" );
1622 : }
1623 :
1624 1 : const realT currentHeight = z_mean();
1625 1 : if( currentHeight == 0 )
1626 : {
1627 0 : if( zm != 0 )
1628 : {
1629 0 : return internal::mxlib_error_report(
1630 0 : error_t::invalidconfig,
1631 0 : "a zero-height atmosphere cannot be rescaled to a positive mean height" );
1632 : }
1633 : }
1634 : else
1635 : {
1636 1 : z_mean( zm );
1637 : }
1638 : }
1639 :
1640 10 : return validate();
1641 20 : }
1642 :
1643 : extern template class aoAtmosphere<float>;
1644 :
1645 : extern template class aoAtmosphere<double>;
1646 :
1647 : extern template class aoAtmosphere<long double>;
1648 :
1649 : #ifdef HASQUAD
1650 : extern template class aoAtmosphere<__float128>;
1651 : #endif
1652 :
1653 : } // namespace analysis
1654 : } // namespace AO
1655 : } // namespace mx
1656 :
1657 : #endif // aoAtmosphere_hpp
|
