Line data Source code
1 : /** \file turbSubHarmonic.hpp
2 : * \brief A class to manage low-frequency sub-harmonic phase screen generation in atmospheric turbulence.
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup mxAO_sim_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2023 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef mx_AO_sim_turbSubHarm_hpp
30 : #define mx_AO_sim_turbSubHarm_hpp
31 :
32 : #include <vector>
33 :
34 : #include "../../base/changeable.hpp"
35 :
36 : #include "../../improc/eigenImage.hpp"
37 : #include "../../improc/eigenCube.hpp"
38 :
39 : #include "../../math/constants.hpp"
40 : #include "../../math/func/jinc.hpp"
41 :
42 : namespace mx
43 : {
44 : namespace AO
45 : {
46 : namespace sim
47 : {
48 :
49 : /// A class to manage low-frequency sub-harmonic phase screen generation in atmospheric turbulence.
50 : /** Implements the method of Johansson and Gavel (1994)\cite johansson_gavel_1994. This is
51 : * needed to generate adequate tip/tilt variance with large outer scales.
52 : *
53 : * \ingroup mxAOSim
54 : */
55 : template <typename _turbAtmosphereT, class _verboseT = mx::verbose::d>
56 : class turbSubHarmonic : public base::changeable<turbSubHarmonic<_turbAtmosphereT>>
57 : {
58 :
59 : public:
60 : typedef _turbAtmosphereT turbAtmosphereT;
61 : typedef _verboseT verboseT;
62 :
63 : typedef typename turbAtmosphereT::realT realT;
64 :
65 : protected:
66 : /** \name Configuration Parameters
67 : * @{
68 : */
69 :
70 : turbAtmosphereT *m_turbAtmo{ nullptr }; ///< Pointer to the parent atmosphere object.
71 :
72 : unsigned m_level{ 1 }; ///< The subharmonic level to apply.
73 :
74 : bool m_outerSubHarmonics{ true }; ///< Whether or not to include the outer subharmonics
75 :
76 : bool m_preCalc{ false }; ///< whether or not the modes are pre-calculated.
77 :
78 : ///@}
79 :
80 : /** \name Internal State
81 : * @{
82 : */
83 :
84 : uint32_t m_scrnSz; ///< The wavefront screen size from the layer being simulated.
85 :
86 : /// Gaussian cosine coefficients prepared for each screen generation.
87 : std::vector<realT> m_cosNoise;
88 :
89 : /// Gaussian sine coefficients prepared for each screen generation.
90 : std::vector<realT> m_sinNoise;
91 :
92 : std::vector<realT> m_m; ///< m-coordinate fractional spatial frequency indices of the subharmonics
93 : std::vector<realT> m_n; ///< n-coordinate fractional spatial frequency indices of the subharmonics
94 :
95 : /// Square root of the PSD-cell variance for each unique conjugate pair.
96 : std::vector<realT> m_sqrtPSD;
97 :
98 : /// Pre-calculated cosine modes followed by their sine counterparts.
99 : improc::eigenCube<realT> m_modes;
100 :
101 : ///@}
102 :
103 : public:
104 : /** \name Construction
105 : * @{
106 : */
107 : turbSubHarmonic();
108 :
109 : ///@}
110 :
111 : /** \name Configuration
112 : * @{
113 : */
114 :
115 : /// Set the pointer to the turbulent atmosphere
116 : void turbAtmo( turbAtmosphereT *atm /**< [in] the new pointer to an AO atmosphere*/ );
117 :
118 : /// Get the pointer to the AO atmosphere
119 : /**
120 : * \returns the the current pointer to the AO atmosphere
121 : */
122 : turbAtmosphereT *turbAtmo();
123 :
124 : /// Set the subharmonic level to apply
125 : void level( uint32_t ml /**< [in] the new level*/ );
126 :
127 : /// Get the subharmonic level
128 : /**
129 : * \returns the current subharmonic level
130 : */
131 : uint32_t level();
132 :
133 : /// Set whether or not the outer subharmonics are included
134 : void outerSubHarmonics( bool osh /**< [in] the new value of the \ref m_outerSubHarmonics flag */ );
135 :
136 : /// Get whether or not the outer subharmonics are included
137 : /**
138 : * \returns the current value of the \ref m_outerSubHarmonics flag
139 : */
140 : bool outerSubHarmonics();
141 :
142 : /// Set whether or not to pre-calculate the modes
143 : void preCalc( bool pc /**< [in] the new value of the \ref m_preCalc flag*/ );
144 :
145 : /// Get whether or not the modes are pre-calculated
146 : /**
147 : * \returns the current value of the preCalc flag
148 : */
149 : bool preCalc();
150 :
151 : ///@}
152 :
153 : /** \name Screen Generation
154 : * @{
155 : */
156 :
157 : /// Allocate needed memory and initialize the subharmonic transform
158 : void initGrid( uint32_t layerNo );
159 :
160 : /// Generate a realization of the subharmonic phase screen and add it to the input screen
161 : void screen( improc::eigenImage<realT> & scrn /**< [in] the input phase screen to which to add the low-frequency screen. Must be m_scrnSz X m_scrnsz*/ );
162 :
163 : /// Deallocate memory
164 : void deinit();
165 :
166 : ///@}
167 : };
168 :
169 : template <typename turbAtmosphereT, class verboseT>
170 4 : turbSubHarmonic<turbAtmosphereT, verboseT>::turbSubHarmonic()
171 : {
172 4 : }
173 :
174 : template <typename turbAtmosphereT, class verboseT>
175 4 : void turbSubHarmonic<turbAtmosphereT, verboseT>::turbAtmo( turbAtmosphereT *turbatm )
176 : {
177 4 : if( turbatm != m_turbAtmo )
178 : {
179 4 : m_turbAtmo = turbatm;
180 4 : this->changed();
181 : }
182 4 : }
183 :
184 : template <typename turbAtmosphereT, class verboseT>
185 : turbAtmosphereT *turbSubHarmonic<turbAtmosphereT, verboseT>::turbAtmo()
186 : {
187 : return m_turbAtmo;
188 : }
189 :
190 : template <typename turbAtmosphereT, class verboseT>
191 4 : void turbSubHarmonic<turbAtmosphereT, verboseT>::level( uint32_t ml )
192 : {
193 4 : if( ml != m_level )
194 : {
195 0 : m_level = ml;
196 0 : this->changed();
197 : }
198 4 : }
199 :
200 : template <typename turbAtmosphereT, class verboseT>
201 : uint32_t turbSubHarmonic<turbAtmosphereT, verboseT>::level()
202 : {
203 : return m_level;
204 : }
205 :
206 : template <typename turbAtmosphereT, class verboseT>
207 4 : void turbSubHarmonic<turbAtmosphereT, verboseT>::outerSubHarmonics( bool osh )
208 : {
209 4 : if( osh != m_outerSubHarmonics )
210 : {
211 0 : m_outerSubHarmonics = osh;
212 0 : this->changed();
213 : }
214 4 : }
215 :
216 : template <typename turbAtmosphereT, class verboseT>
217 : bool turbSubHarmonic<turbAtmosphereT, verboseT>::outerSubHarmonics()
218 : {
219 : return m_outerSubHarmonics;
220 : }
221 :
222 : template <typename turbAtmosphereT, class verboseT>
223 4 : void turbSubHarmonic<turbAtmosphereT, verboseT>::preCalc( bool pc )
224 : {
225 4 : if( pc != m_preCalc )
226 : {
227 3 : m_preCalc = pc;
228 3 : this->changed();
229 : }
230 4 : }
231 :
232 : template <typename turbAtmosphereT, class verboseT>
233 : bool turbSubHarmonic<turbAtmosphereT, verboseT>::preCalc()
234 : {
235 : return m_preCalc;
236 : }
237 :
238 : template <typename turbAtmosphereT, class verboseT>
239 4 : void turbSubHarmonic<turbAtmosphereT, verboseT>::initGrid( uint32_t layerNo )
240 : {
241 4 : if( m_turbAtmo == nullptr )
242 : {
243 0 : throw mx::exception<verboseT>( error_t::paramnotset, "atmosphere is not set (m_turbAtmo is nullptr)" );
244 : }
245 :
246 4 : if( m_turbAtmo->aosys() == nullptr )
247 : {
248 0 : throw mx::exception<verboseT>( error_t::paramnotset, "ao system is not set (m_turbAtmo->m_aosys is nullptr)" );
249 : }
250 :
251 4 : if( m_turbAtmo->nLayers() <= layerNo )
252 : {
253 0 : throw mx::exception<verboseT>( error_t::invalidconfig,
254 0 : "atmosphere is not setup (m_turbAtmo->m_layers size is <= layerNo)" );
255 : }
256 :
257 4 : m_m.resize( 0 ); // resized dynamically below
258 4 : m_n.resize( 0 );
259 4 : m_sqrtPSD.resize( 0 );
260 4 : m_cosNoise.resize( 0 );
261 4 : m_sinNoise.resize( 0 );
262 :
263 4 : m_scrnSz = m_turbAtmo->layer( layerNo ).scrnSz();
264 :
265 4 : if( m_level == 0 )
266 : {
267 :
268 0 : m_modes.clear();
269 0 : this->setChangePoint();
270 0 : return;
271 : }
272 :
273 4 : realT r0 = m_turbAtmo->aosys()->atm.r_0( m_turbAtmo->aosys()->lam_sci() );
274 4 : realT D = m_turbAtmo->aosys()->D();
275 4 : uint32_t wfSz = m_turbAtmo->wfSz();
276 :
277 4 : realT beta = 0.0218 / pow( r0, 5. / 3. ) / pow( ( D / wfSz ) * m_scrnSz, 2 );
278 :
279 4 : realT sqrt_alpha = 0.5 * 11. / 3.;
280 :
281 4 : realT L0 = m_turbAtmo->aosys()->atm.L_0( layerNo );
282 :
283 : realT L02;
284 4 : if( L0 > 0 )
285 0 : L02 = 1.0 / pow( L0, 2 );
286 : else
287 4 : L02 = 0;
288 :
289 4 : std::vector<realT> scs;
290 :
291 4 : if( m_outerSubHarmonics )
292 : {
293 4 : realT sc = 0.25;
294 :
295 8 : m_m = std::vector<realT>( { -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, -1.25, 1.25, -1.25, 1.25,
296 : -1.25, 1.25, -1.25, 1.25, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25 } );
297 4 : m_n = std::vector<realT>( { -1.25, -1.25, -1.25, -1.25, -1.25, -1.25, -0.75, -0.75, -0.25, -0.25,
298 : 0.25, 0.25, 0.75, 0.75, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25 } );
299 4 : scs.resize( m_m.size(), sc );
300 : }
301 :
302 8 : for( int nl = 1; nl <= m_level; ++nl )
303 : {
304 4 : realT sc = pow( 3.0, -nl );
305 :
306 28 : for( int mp = -3; mp < 3; ++mp )
307 : {
308 168 : for( int np = -3; np < 3; ++np )
309 : {
310 144 : if( nl < m_level )
311 : {
312 0 : if( mp == -1 )
313 : {
314 0 : if( np == -1 || np == 0 )
315 0 : continue;
316 : }
317 0 : else if( mp == 0 )
318 : {
319 0 : if( np == -1 || np == 0 )
320 0 : continue;
321 : }
322 : }
323 :
324 144 : m_m.push_back( sc * ( mp + 0.5 ) );
325 144 : m_n.push_back( sc * ( np + 0.5 ) );
326 144 : scs.push_back( sc );
327 : }
328 : }
329 : }
330 :
331 : // Retain one member of each +/-k pair. A real phase screen requires both
332 : // quadratures for that pair, rather than two identical cosine modes.
333 4 : std::vector<realT> uniqueM;
334 4 : std::vector<realT> uniqueN;
335 4 : std::vector<realT> uniqueScs;
336 228 : for( size_t n = 0; n < m_m.size(); ++n )
337 : {
338 224 : if( m_m[n] > 0 || ( m_m[n] == 0 && m_n[n] > 0 ) )
339 : {
340 112 : uniqueM.push_back( m_m[n] );
341 112 : uniqueN.push_back( m_n[n] );
342 112 : uniqueScs.push_back( scs[n] );
343 : }
344 : }
345 :
346 4 : m_m.swap( uniqueM );
347 4 : m_n.swap( uniqueN );
348 4 : scs.swap( uniqueScs );
349 4 : m_sqrtPSD.resize( m_m.size() );
350 4 : m_cosNoise.resize( m_m.size() );
351 4 : m_sinNoise.resize( m_m.size() );
352 :
353 4 : if( m_preCalc )
354 : {
355 3 : m_modes.resize( m_scrnSz, m_scrnSz, 2 * m_m.size() );
356 : }
357 :
358 116 : for( size_t n = 0; n < m_m.size(); ++n )
359 : {
360 112 : realT k = sqrt( ( pow( m_m[n], 2 ) + pow( m_n[n], 2 ) ) ) / ( ( D / wfSz ) * m_scrnSz );
361 :
362 112 : realT tpsd = beta / pow( k * k + L02, sqrt_alpha );
363 :
364 112 : realT Ppiston = 0;
365 112 : realT Ptiptilt = 0;
366 112 : if( m_turbAtmo->aosys()->psd.subPiston() )
367 : {
368 112 : Ppiston = pow( 2 * math::func::jinc( math::pi<realT>() * k * D ), 2 );
369 : }
370 :
371 112 : if( m_turbAtmo->aosys()->psd.subTipTilt() )
372 : {
373 0 : Ptiptilt = pow( 4 * math::func::jincN( 2, math::pi<realT>() * k * D ), 2 );
374 : }
375 :
376 : // The representative stands for both members of its conjugate pair.
377 : // The factor sqrt(2) gives the pair's total PSD-cell variance when
378 : // independent unit-variance cosine and sine coefficients are used.
379 112 : m_sqrtPSD[n] = sqrt( 2 ) * scs[n] * sqrt( tpsd * ( 1 - Ppiston - Ptiptilt ) );
380 :
381 112 : if( m_preCalc )
382 : {
383 9044 : for( int cc = 0; cc < m_scrnSz; ++cc )
384 : {
385 8960 : realT np = cc - 0.5 * m_scrnSz;
386 1270528 : for( int rr = 0; rr < m_scrnSz; ++rr )
387 : {
388 1261568 : realT mp = rr - 0.5 * m_scrnSz;
389 1261568 : m_modes.image( n )( rr, cc ) =
390 1261568 : m_sqrtPSD[n] * cos( math::two_pi<realT>() * ( m_m[n] * mp + m_n[n] * np ) / m_scrnSz );
391 1261568 : m_modes.image( n + m_m.size() )( rr, cc ) =
392 1261568 : m_sqrtPSD[n] * sin( math::two_pi<realT>() * ( m_m[n] * mp + m_n[n] * np ) / m_scrnSz );
393 : }
394 : }
395 : }
396 : }
397 :
398 4 : this->setChangePoint();
399 4 : }
400 :
401 : template <typename turbAtmosphereT, class verboseT>
402 11 : void turbSubHarmonic<turbAtmosphereT, verboseT>::screen( improc::eigenImage<realT> &scrn )
403 : {
404 11 : if( m_level == 0 )
405 : {
406 0 : return;
407 : }
408 :
409 11 : if( m_turbAtmo == nullptr )
410 : {
411 0 : throw mx::exception<verboseT>( error_t::paramnotset, "atmosphere is not set (m_turbAtmo is nullptr)" );
412 : }
413 :
414 11 : if( this->isChanged() )
415 : {
416 0 : throw mx::exception<verboseT>( error_t::invalidconfig, "configuration has changed but not re-initialized" );
417 : }
418 :
419 11 : if( scrn.rows() != m_scrnSz || scrn.cols() != m_scrnSz )
420 : {
421 0 : throw mx::exception<verboseT>( error_t::sizeerr, "input screen is not the right size" );
422 : }
423 :
424 : // Check that we're allocated
425 11 : if( m_preCalc )
426 : {
427 20 : if( m_modes.rows() != scrn.rows() || m_modes.cols() != scrn.cols() ||
428 10 : m_modes.planes() != 2 * m_cosNoise.size() )
429 : {
430 0 : throw mx::exception<verboseT>( error_t::sizeerr, "modes cube wrong size, call initGrid()." );
431 : }
432 : }
433 : else
434 : {
435 2 : if( m_cosNoise.size() != m_m.size() || m_sinNoise.size() != m_m.size() || m_m.size() != m_n.size() ||
436 1 : m_sqrtPSD.size() != m_m.size() )
437 : {
438 0 : throw mx::exception<verboseT>( error_t::sizeerr, "vectors not allocated, call initGrid()." );
439 : }
440 : }
441 :
442 : // Now fill in the noise
443 319 : for( size_t n = 0; n < m_cosNoise.size(); ++n )
444 : {
445 308 : m_cosNoise[n] = m_turbAtmo->normVar();
446 308 : m_sinNoise[n] = m_turbAtmo->normVar();
447 : }
448 :
449 11 : #pragma omp parallel for
450 : for( int cc = 0; cc < m_scrnSz; ++cc )
451 : {
452 : realT np = cc - 0.5 * m_scrnSz;
453 : for( int rr = 0; rr < m_scrnSz; ++rr )
454 : {
455 : realT mp = rr - 0.5 * m_scrnSz;
456 :
457 : if( m_preCalc )
458 : {
459 : for( unsigned n = 0; n < m_cosNoise.size(); ++n )
460 : {
461 : scrn( rr, cc ) += m_cosNoise[n] * m_modes.image( n )( rr, cc ) +
462 : m_sinNoise[n] * m_modes.image( n + m_cosNoise.size() )( rr, cc );
463 : }
464 : }
465 : else
466 : {
467 : for( unsigned n = 0; n < m_m.size(); ++n )
468 : {
469 : realT arg = math::two_pi<realT>() * ( m_m[n] * mp + m_n[n] * np ) / m_scrnSz;
470 : scrn( rr, cc ) += m_sqrtPSD[n] * ( m_cosNoise[n] * cos( arg ) + m_sinNoise[n] * sin( arg ) );
471 : }
472 : }
473 : }
474 : }
475 : }
476 :
477 : template <typename turbAtmosphereT, class verboseT>
478 0 : void turbSubHarmonic<turbAtmosphereT, verboseT>::deinit()
479 : {
480 0 : m_m.clear();
481 0 : m_n.clear();
482 0 : m_sqrtPSD.clear();
483 0 : m_cosNoise.clear();
484 0 : m_sinNoise.clear();
485 0 : m_modes.resize( 0, 0, 0 );
486 :
487 0 : this->changed();
488 0 : }
489 :
490 : } // namespace sim
491 : } // namespace AO
492 : } // namespace mx
493 :
494 : #endif // mx_AO_sim_turbSubHarm_hpp
|
