Line data Source code
1 : /** \file turbAtmosphere.hpp
2 : * \brief Declaration and definition of a turbulent atmosphere.
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_turbAtmosphere_hpp
30 : #define mx_AO_sim_turbAtmosphere_hpp
31 :
32 : #include <vector>
33 : #include <iostream>
34 :
35 : #include "../../sigproc/psdFilter.hpp"
36 : #include "../../sigproc/psdUtils.hpp"
37 :
38 : #include "../../base/changeable.hpp"
39 :
40 : #include "../../improc/milkImage.hpp"
41 :
42 : #include "../../math/constants.hpp"
43 : #include "../../math/func/jinc.hpp"
44 : #include "../../math/randomT.hpp"
45 :
46 : #include "../../ioutils/fits/fitsFile.hpp"
47 : #include "../../ioutils/stringUtils.hpp"
48 : #include "../../sys/timeUtils.hpp"
49 :
50 : #include "turbLayer.hpp"
51 : #include "turbSubHarmonic.hpp"
52 : #include "wavefront.hpp"
53 :
54 : #ifdef DEBUG
55 : #define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
56 : #else
57 : #define BREAD_CRUMB
58 : #endif
59 :
60 : namespace mx
61 : {
62 : namespace AO
63 : {
64 : namespace sim
65 : {
66 :
67 : /// A turbulent atmosphere simulator
68 : /** Generate multi-layer phase screens using PSD filtering with the FFT. Manages
69 : * the shifting and combination of the phase screens (without propagation).
70 : *
71 : * The \ref turbSubHarmonic class provides low-frequency sub-harmonics if desired
72 : *
73 : * The \ref turbLayer class manages the actual interpolation.
74 : *
75 : * \todo add facility for switching interpolators
76 : * \todo need to include wrap detection when using subharmonics (e.g. frames needs to be right)
77 : *
78 : * \ingroup mxAOSim
79 : */
80 : template <typename _aoSystemT, class _verboseT = mx::verbose::d>
81 : struct turbAtmosphere : public base::changeable<turbAtmosphere<_aoSystemT, _verboseT>>
82 : {
83 :
84 : public:
85 : typedef _aoSystemT aoSystemT;
86 : typedef _verboseT verboseT;
87 :
88 : typedef typename aoSystemT::realT realT;
89 : typedef Eigen::Array<realT, -1, -1> imageT;
90 :
91 : protected:
92 : /** \name Configuration Data
93 : * @{
94 : */
95 : uint32_t m_wfSz{ 0 }; ///< Size of the wavefront in pixels.
96 :
97 : uint32_t m_buffSz{ 0 }; ///< Buffer to apply around wavefront for interpolation.
98 :
99 : aoSystemT *m_aosys{ nullptr }; ///< The AO system object containing all system parameters.
100 :
101 : uint32_t m_shLevel{ 0 }; /**< Number of subharmonic levels. 0 means no subharmonics. Generally only
102 : * 1 level is needed to obtain good Tip/Tilt results. See \ref turbSubHarmonic.
103 : */
104 :
105 : bool m_outerSubHarmonics{ true }; /**< Whether or not the outer subharmonics are included.
106 : * See \ref turbSubHarmonic.
107 : */
108 :
109 : bool m_shPreCalc{ true }; /**< Whether or not to pre-calculate the subharmonic modes. This is a
110 : * trade of speed vs. memory use. See \ref turbSubHarmonic.
111 : */
112 :
113 : bool m_retain{ false }; /**< Whether or not to retain working memory after screen generation. One
114 : * would set to true if many screens will be calculated in monte carlo fashion.
115 : * For a standard case of one set of screens being shifted by wind velocity,
116 : * this should be false to minimize memory use.
117 : */
118 :
119 : bool m_forceGen{ true }; /**< Force generation of new screens if true. Note that this class does not presently
120 : * do any checks of saved screens to verify they match the current configuration.
121 : * Set to true with caution.
122 : *
123 : * \todo need check for changes in parameters to decide if saved frames are valid per layer
124 : */
125 :
126 : std::string m_dataDir; /**< Specifies a path where phase screens are stored. Leaving this unset or "" is equivalent
127 : * to \ref m_forceGen` == true`.
128 : */
129 :
130 : imageT *m_pupil{ 0 }; ///< A pointer to the pupil mask.
131 :
132 : realT m_timeStep{ 0 }; ///< Length of each iteration, in seconds.
133 :
134 : ///@}
135 :
136 : /** \name Internal State
137 : * @{
138 : */
139 : std::vector<turbLayer<aoSystemT>> m_layers; ///< Vector of turbulent layers.
140 :
141 : std::vector<turbSubHarmonic<turbAtmosphere>> m_subHarms; ///< Sub-harmonic layers
142 :
143 : size_t m_frames{ 0 }; ///< Length of the turbulence sequence.
144 :
145 : int m_nWf{ 0 }; ///< Number of iterations which have occurred.
146 :
147 : math::normDistT<realT> m_normVar; ///< Normal random deviate generator. This seeded in the constructor.
148 :
149 : public:
150 : /** \name Construction and Configuration
151 : * @{
152 : */
153 :
154 : /// Default c'tor
155 : turbAtmosphere();
156 :
157 : /// Setup the overall atmosphere.
158 : /**
159 : */
160 : void setup( uint32_t wfSz, ///< [in] The size of the wavefront in pixels.
161 : uint32_t buffSz, ///< [in] The size of the interpolation buffer to use.
162 : aoSystemT *aosys, ///< [in] Pointer to an AO System. See \ref m_aosys.
163 : uint32_t shLevel /**< [in] number of subharmonic levels to use. 0 turns off subharmonics.
164 : See \ref m_shLevel.*/
165 : );
166 :
167 : /// Set the wavefront size
168 : /**
169 : * see \ref m_wfSz
170 : */
171 : void wfSz( uint32_t ws /**< [in] the new wavefront size */ );
172 :
173 : /// Get the wavefront size
174 : /**
175 : * \returns the current value of \ref m_wfSz
176 : */
177 : uint32_t wfSz();
178 :
179 : /// Set the interpolation buffer size
180 : /**
181 : * see \ref m_buffSz
182 : */
183 : void buffSz( uint32_t bs /**< [in] the new buffer size*/ );
184 :
185 : /// Get the interpolation buffer size
186 : /**
187 : * \returns the current value of \ref m_buffSz
188 : */
189 : uint32_t buffSz();
190 :
191 : /// Set the pointer to an AO system object
192 : /**
193 : * see \ref m_aosys
194 : */
195 : void aosys( aoSystemT *aos /**< [in] the new pointer to an AO system object*/ );
196 :
197 : /// Get the pointer to an AO system object
198 : /**
199 : * \returns the current value of \ref m_aosys
200 : */
201 : aoSystemT *aosys();
202 :
203 : /// Set the subharmonic level
204 : /**
205 : * see \ref m_shLevel
206 : */
207 : void shLevel( uint32_t shl /**< [in] the new subharmonic level*/ );
208 :
209 : /// Get the subharmonic level
210 : /**
211 : * \returns the current value of \ref m_shLevel
212 : */
213 : uint32_t shLevel();
214 :
215 : /// Set whether outer subharmonics are included
216 : /**
217 : * see \ref m_outerSubHarmonics
218 : */
219 : void outerSubHarmonics( bool osh /**< [in] the new value of flag controlling outer subharmonics */ );
220 :
221 : /// Get whether outer subharmonics are included
222 : /**
223 : * \returns the current value of \ref m_outerSubHarmonics
224 : */
225 : bool outerSubHarmonics();
226 :
227 : /// Set whether subharmonic modes are pre-calculated
228 : /**
229 : * see \ref m_shPreCalc
230 : */
231 : void shPreCalc( bool shp /**< [in] the new value of flag controlling subharmonic mode precalculation */ );
232 :
233 : /// Get whether subharmonic modes are pre-calculated
234 : /**
235 : * \returns the current value of \ref m_shPreCalc
236 : */
237 : bool shPreCalc();
238 :
239 : /// Set whether memory for screen generation is retained
240 : /**
241 : * see \ref m_retain
242 : */
243 : void retain( bool rtn /**< [in] the new value of the flag controlling whether memory is retained*/ );
244 :
245 : /// Get whether memory for screen generation is retained
246 : /**
247 : * \returns the current value of \ref m_retain
248 : */
249 : bool retain();
250 :
251 : /// Set whether new screen generation is forced
252 : /**
253 : * see \ref m_forceGen
254 : */
255 : void forceGen( bool fg /**< [in] the new value of the flag controlling whether screen generation is forced*/ );
256 :
257 : /// Get whether new screen generation is forced
258 : /**
259 : * \returns the current value of m_forceGen
260 : */
261 : bool forceGen();
262 :
263 : /// Set the data directory for saving phase screens
264 : /**
265 : * see \ref m_dataDir
266 : */
267 : void dataDir( const std::string &dd /**< [in] the new data directory for phase screen saving*/ );
268 :
269 : /// Get the data directory for saving phase screens
270 : /**
271 : * \returns the current value of m_dataDir
272 : */
273 : std::string dataDir();
274 :
275 : /// Setup the layers and prepare them for phase screen generation.
276 : /** The number of entries in \p scrnSz must mach the number of layers in the atmosphere
277 : * of \ref m_aosys.
278 : */
279 : void setLayers( const std::vector<size_t> &scrnSz /**< [in] the screen sizes for each layer.*/ );
280 :
281 : /// Setup the layers and prepare them for phase screen generation.
282 : /** This sets all layer phase screens to be the same size.
283 : *
284 : */
285 : void setLayers( const size_t scrnSz /**< [in] the screen sizes for all layers.*/ );
286 :
287 : /// Get the number of layers
288 : /**
289 : * \returns the size of the layers vector \ref m_layers
290 : */
291 : uint32_t nLayers();
292 :
293 : /// Get a layer
294 : /**
295 : * \returns a reference to on entry in \ref m_layers
296 : */
297 : turbLayer<aoSystemT> &layer( uint32_t n );
298 :
299 : /// Get the random number generator
300 : /**
301 : * \returns a reference to \ref m_normVar
302 : */
303 : math::normDistT<realT> &normVar();
304 :
305 : ///@} - Construction and Configuration
306 :
307 : uint32_t scrnLengthPixels();
308 :
309 : uint32_t scrnLengthPixels( uint32_t n );
310 :
311 : realT scrnLength();
312 :
313 : realT scrnLength( uint32_t n );
314 :
315 : uint32_t maxShift( realT dt );
316 :
317 : uint32_t maxShift( uint32_t n, realT dt );
318 :
319 : /** \name Screen Generation
320 : * @{
321 : */
322 :
323 : /// Generate all phase screens
324 : /** Loads them from disk if possible and \ref m_forceGen == false.
325 : *
326 : * Deallocates working memory when done, unless \ref m_retain == true.
327 : */
328 : void genLayers();
329 :
330 : /// @}
331 :
332 : int shift( improc::milkImage<realT> &phase, realT dt );
333 :
334 : int frames( int f );
335 : size_t frames();
336 :
337 : int wfPS( realT ps ); ///< dummy function for simulatedAOSystem. Does nothing.
338 :
339 : bool _loopClosed;
340 :
341 : void nextWF( wavefront<realT> &wf );
342 : };
343 :
344 : template <typename aoSystemT, class verboseT>
345 2 : turbAtmosphere<aoSystemT, verboseT>::turbAtmosphere()
346 : {
347 2 : m_normVar.seed();
348 2 : }
349 :
350 : template <typename aoSystemT, class verboseT>
351 2 : void turbAtmosphere<aoSystemT, verboseT>::setup( uint32_t ws, uint32_t bs, aoSystemT *aos, uint32_t shl )
352 : {
353 2 : wfSz( ws );
354 2 : buffSz( bs );
355 2 : aosys( aos );
356 2 : shLevel( shl );
357 2 : }
358 :
359 : template <typename aoSystemT, class verboseT>
360 2 : void turbAtmosphere<aoSystemT, verboseT>::wfSz( uint32_t ws )
361 : {
362 2 : if( ws != m_wfSz )
363 : {
364 2 : m_wfSz = ws;
365 2 : this->changed();
366 : }
367 2 : }
368 :
369 : template <typename aoSystemT, class verboseT>
370 8 : uint32_t turbAtmosphere<aoSystemT, verboseT>::wfSz()
371 : {
372 8 : return m_wfSz;
373 : }
374 :
375 : template <typename aoSystemT, class verboseT>
376 2 : void turbAtmosphere<aoSystemT, verboseT>::buffSz( uint32_t bs )
377 : {
378 2 : if( bs != m_buffSz )
379 : {
380 0 : m_buffSz = bs;
381 0 : this->changed();
382 : }
383 2 : }
384 :
385 : template <typename aoSystemT, class verboseT>
386 2 : uint32_t turbAtmosphere<aoSystemT, verboseT>::buffSz()
387 : {
388 2 : return m_buffSz;
389 : }
390 :
391 : template <typename aoSystemT, class verboseT>
392 2 : void turbAtmosphere<aoSystemT, verboseT>::aosys( aoSystemT *aos )
393 : {
394 2 : if( aos != m_aosys )
395 : {
396 2 : m_aosys = aos;
397 2 : this->changed();
398 : }
399 2 : }
400 :
401 : template <typename aoSystemT, class verboseT>
402 123144 : aoSystemT *turbAtmosphere<aoSystemT, verboseT>::aosys()
403 : {
404 123144 : return m_aosys;
405 : }
406 :
407 : template <typename aoSystemT, class verboseT>
408 2 : void turbAtmosphere<aoSystemT, verboseT>::shLevel( uint32_t shl )
409 : {
410 2 : if( shl != m_shLevel )
411 : {
412 2 : m_shLevel = shl;
413 2 : this->changed();
414 : }
415 2 : }
416 :
417 : template <typename aoSystemT, class verboseT>
418 2 : uint32_t turbAtmosphere<aoSystemT, verboseT>::shLevel()
419 : {
420 2 : return m_shLevel;
421 : }
422 :
423 : template <typename aoSystemT, class verboseT>
424 : void turbAtmosphere<aoSystemT, verboseT>::outerSubHarmonics( bool osh )
425 : {
426 : if( osh != m_outerSubHarmonics )
427 : {
428 : m_outerSubHarmonics = osh;
429 : this->changed();
430 : }
431 : }
432 :
433 : template <typename aoSystemT, class verboseT>
434 2 : bool turbAtmosphere<aoSystemT, verboseT>::outerSubHarmonics()
435 : {
436 2 : return m_outerSubHarmonics;
437 : }
438 :
439 : template <typename aoSystemT, class verboseT>
440 : void turbAtmosphere<aoSystemT, verboseT>::shPreCalc( bool shp )
441 : {
442 : if( shp != m_shPreCalc )
443 : {
444 : m_shPreCalc = shp;
445 : this->changed();
446 : }
447 : }
448 :
449 : template <typename aoSystemT, class verboseT>
450 : bool turbAtmosphere<aoSystemT, verboseT>::shPreCalc()
451 : {
452 : return m_shPreCalc;
453 : }
454 :
455 : template <typename aoSystemT, class verboseT>
456 1 : void turbAtmosphere<aoSystemT, verboseT>::retain( bool rtn )
457 : {
458 1 : if( rtn != m_retain )
459 : {
460 1 : m_retain = rtn;
461 1 : this->changed();
462 : }
463 1 : }
464 :
465 : template <typename aoSystemT, class verboseT>
466 : bool turbAtmosphere<aoSystemT, verboseT>::retain()
467 : {
468 : return m_retain;
469 : }
470 :
471 : template <typename aoSystemT, class verboseT>
472 : void turbAtmosphere<aoSystemT, verboseT>::forceGen( bool fg )
473 : {
474 : if( fg != m_forceGen )
475 : {
476 : this->changed();
477 : }
478 :
479 : m_forceGen = fg;
480 : }
481 :
482 : template <typename aoSystemT, class verboseT>
483 : bool turbAtmosphere<aoSystemT, verboseT>::forceGen()
484 : {
485 : return m_forceGen;
486 : }
487 :
488 : template <typename aoSystemT, class verboseT>
489 : void turbAtmosphere<aoSystemT, verboseT>::dataDir( const std::string &dd )
490 : {
491 : if( dd != m_dataDir )
492 : {
493 : this->changed();
494 : }
495 :
496 : m_dataDir = dd;
497 : }
498 :
499 : template <typename aoSystemT, class verboseT>
500 : std::string turbAtmosphere<aoSystemT, verboseT>::dataDir()
501 : {
502 : return m_dataDir;
503 : }
504 :
505 : template <typename aoSystemT, class verboseT>
506 2 : void turbAtmosphere<aoSystemT, verboseT>::setLayers( const std::vector<size_t> &scrnSz )
507 : {
508 2 : if( m_aosys == nullptr )
509 : {
510 0 : throw mx::exception( error_t::paramnotset, "ao system is not set (m_aosys is nullptr)" );
511 : }
512 :
513 2 : size_t nLayers = scrnSz.size();
514 :
515 2 : if( nLayers != m_aosys->atm.n_layers() )
516 : {
517 0 : throw mx::exception( error_t::invalidarg,
518 :
519 0 : "Size of scrnSz vector does not match atmosphere." );
520 : }
521 :
522 2 : if( nLayers != m_layers.size() )
523 : {
524 2 : this->changed();
525 : }
526 :
527 2 : m_layers.resize( nLayers );
528 :
529 4 : for( size_t i = 0; i < nLayers; ++i )
530 : {
531 2 : m_layers[i].setLayer( this, i, scrnSz[i] );
532 : }
533 :
534 2 : if( m_shLevel > 0 )
535 : {
536 2 : if( nLayers != m_subHarms.size() )
537 : {
538 2 : this->changed();
539 : }
540 :
541 2 : m_subHarms.resize( nLayers );
542 :
543 4 : for( size_t i = 0; i < nLayers; ++i )
544 : {
545 2 : m_subHarms[i].turbAtmo( this );
546 2 : m_subHarms[i].preCalc( m_shPreCalc );
547 2 : m_subHarms[i].level( m_shLevel );
548 2 : m_subHarms[i].outerSubHarmonics( m_outerSubHarmonics );
549 2 : m_subHarms[i].initGrid( i );
550 : }
551 : }
552 2 : }
553 :
554 : template <typename aoSystemT, class verboseT>
555 2 : void turbAtmosphere<aoSystemT, verboseT>::setLayers( const size_t scrnSz )
556 : {
557 2 : if( m_aosys == nullptr )
558 : {
559 0 : throw mx::exception( error_t::paramnotset,
560 :
561 0 : "atmosphere is not set (m_atm is nullptr)" );
562 : }
563 :
564 2 : size_t n = m_aosys->atm.n_layers();
565 :
566 2 : setLayers( std::vector<size_t>( n, scrnSz ) );
567 2 : }
568 :
569 : template <typename aoSystemT, class verboseT>
570 5 : uint32_t turbAtmosphere<aoSystemT, verboseT>::nLayers()
571 : {
572 5 : return m_layers.size();
573 : }
574 :
575 : template <typename aoSystemT, class verboseT>
576 6 : turbLayer<aoSystemT> &turbAtmosphere<aoSystemT, verboseT>::layer( uint32_t n )
577 : {
578 6 : if( n >= m_layers.size() )
579 : {
580 0 : throw mx::exception( error_t::invalidarg, "n too large for number of layers." );
581 : }
582 :
583 6 : return m_layers[n];
584 : }
585 :
586 : template <typename aoSystemT, class verboseT>
587 619 : math::normDistT<typename aoSystemT::realT> &turbAtmosphere<aoSystemT, verboseT>::normVar()
588 : {
589 619 : return m_normVar;
590 : }
591 :
592 : template <typename aoSystemT, class verboseT>
593 : uint32_t turbAtmosphere<aoSystemT, verboseT>::scrnLengthPixels()
594 : {
595 : return 0;
596 : }
597 :
598 : template <typename aoSystemT, class verboseT>
599 : uint32_t turbAtmosphere<aoSystemT, verboseT>::scrnLengthPixels( uint32_t n )
600 : {
601 :
602 : return 0;
603 : /* realT dx = (scrnSz-wfSz-turb.buffSz()) / cos(m_aosys->atm.layer_v_wind(n));
604 :
605 : pow(2*pow(scrnSz-wfSz-turb.buffSz(),2), 0.5) - 1;
606 : */
607 : }
608 :
609 : template <typename aoSystemT, class verboseT>
610 : turbAtmosphere<aoSystemT, verboseT>::realT turbAtmosphere<aoSystemT, verboseT>::scrnLength()
611 : {
612 : return 0;
613 : }
614 :
615 : template <typename aoSystemT, class verboseT>
616 : turbAtmosphere<aoSystemT, verboseT>::realT turbAtmosphere<aoSystemT, verboseT>::scrnLength( uint32_t n )
617 : {
618 : return 0;
619 : }
620 :
621 : template <typename aoSystemT, class verboseT>
622 : uint32_t turbAtmosphere<aoSystemT, verboseT>::maxShift( realT dt )
623 : {
624 : return 0;
625 : }
626 :
627 : template <typename aoSystemT, class verboseT>
628 : uint32_t turbAtmosphere<aoSystemT, verboseT>::maxShift( uint32_t n, realT dt )
629 : {
630 : return 0;
631 : }
632 :
633 : template <typename aoSystemT, class verboseT>
634 1 : void turbAtmosphere<aoSystemT, verboseT>::genLayers()
635 : {
636 1 : if( m_dataDir != "" && !m_forceGen )
637 : {
638 0 : std::string fbase = m_dataDir;
639 0 : fbase += "/";
640 0 : fbase += "layer_";
641 :
642 0 : fits::fitsFile<realT> ff;
643 0 : std::string fname;
644 0 : for( size_t i = 0; i < m_layers.size(); ++i )
645 : {
646 0 : fname = fbase + ioutils::convertToString<int>( i ) + ".fits";
647 0 : ff.read( m_layers[i].m_phase, fname );
648 : }
649 :
650 0 : this->setChangePoint();
651 :
652 0 : return;
653 0 : }
654 :
655 1 : sigproc::psdFilter<realT, 2> filt;
656 :
657 2 : for( size_t i = 0; i < m_layers.size(); ++i )
658 : {
659 1 : if( m_layers[i].m_phase.rows() != m_layers[i].scrnSz() || m_layers[i].m_phase.cols() != m_layers[i].scrnSz() )
660 : {
661 0 : throw mx::exception( error_t::sizeerr, "layer phase not allocated." );
662 : }
663 :
664 1 : if( m_layers[i].m_freq.rows() != m_layers[i].scrnSz() || m_layers[i].m_freq.cols() != m_layers[i].scrnSz() )
665 : {
666 0 : throw mx::exception( error_t::sizeerr, "layer freq not allocated." );
667 : }
668 :
669 1 : if( m_layers[i].m_psd.rows() != m_layers[i].scrnSz() || m_layers[i].m_psd.cols() != m_layers[i].scrnSz() )
670 : {
671 0 : throw mx::exception( error_t::sizeerr, "layer psd not allocated." );
672 : }
673 :
674 193 : for( size_t jj = 0; jj < m_layers[i].m_scrnSz; ++jj )
675 : {
676 37056 : for( size_t ii = 0; ii < m_layers[i].m_scrnSz; ++ii )
677 : {
678 36864 : m_layers[i].m_phase( ii, jj ) = m_normVar;
679 : }
680 : }
681 :
682 1 : realT dkx = m_layers[i].m_freq( 1, 0 ) - m_layers[i].m_freq( 0, 0 );
683 1 : realT dky = m_layers[i].m_freq( 0, 1 ) - m_layers[i].m_freq( 0, 0 );
684 1 : filt.psdSqrt( m_layers[i].m_psd, dkx, dky );
685 :
686 1 : filt( m_layers[i].m_phase );
687 :
688 1 : if( m_shLevel > 0 )
689 : {
690 1 : m_subHarms[i].screen( m_layers[i].m_phase );
691 : }
692 : }
693 :
694 1 : if( !m_retain )
695 : {
696 0 : for( size_t i = 0; i < m_layers.size(); ++i )
697 : {
698 0 : m_layers[i].genDealloc();
699 : }
700 :
701 0 : for( size_t i = 0; i < m_subHarms.size(); ++i )
702 : {
703 0 : m_subHarms[i].deinit();
704 : }
705 : }
706 :
707 1 : if( m_dataDir != "" )
708 : {
709 0 : std::string fbase = m_dataDir;
710 0 : fbase += "/";
711 0 : fbase += "layer_";
712 :
713 0 : fits::fitsFile<realT> ff;
714 0 : std::string fname;
715 0 : for( size_t i = 0; i < m_layers.size(); ++i )
716 : {
717 0 : fname = fbase + ioutils::convertToString<int>( i ) + ".fits";
718 0 : ff.write( fname, m_layers[i].m_phase );
719 : }
720 0 : }
721 :
722 1 : this->setChangePoint();
723 1 : }
724 :
725 : template <typename aoSystemT, class verboseT>
726 : int turbAtmosphere<aoSystemT, verboseT>::shift( improc::milkImage<realT> &milkPhase, realT dt )
727 : {
728 : if( this->isChanged() )
729 : {
730 : throw mx::exception( error_t::invalidconfig, "configuration has changed but genLayers has not been run." );
731 : }
732 :
733 : improc::eigenMap<realT> phase( milkPhase );
734 :
735 : if( phase.rows() != m_wfSz || phase.cols() != m_wfSz )
736 : {
737 : }
738 :
739 : // Don't use OMP if not multiple layers b/c it seems to make it use multiple threads in some awful way
740 : if( m_layers.size() > 1 )
741 : {
742 : // clang-format off
743 : #pragma omp parallel for // clang-format on
744 : for( size_t j = 0; j < m_layers.size(); ++j )
745 : {
746 : m_layers[j].shift( dt );
747 : }
748 : }
749 : else
750 : {
751 : m_layers[0].shift( dt );
752 : }
753 :
754 : milkPhase.setWrite();
755 : phase.setZero();
756 : for( size_t j = 0; j < m_layers.size(); ++j )
757 : {
758 : phase +=
759 : sqrt( m_aosys->atm.layer_Cn2( j ) ) * m_layers[j].m_shiftPhase.block( m_buffSz, m_buffSz, m_wfSz, m_wfSz );
760 : }
761 : milkPhase.post();
762 :
763 : return 0;
764 : }
765 :
766 : template <typename aoSystemT, class verboseT>
767 : int turbAtmosphere<aoSystemT, verboseT>::frames( int f )
768 : {
769 : m_frames = f;
770 :
771 : return 0;
772 : }
773 :
774 : template <typename aoSystemT, class verboseT>
775 : size_t turbAtmosphere<aoSystemT, verboseT>::frames()
776 : {
777 : return m_frames;
778 : }
779 :
780 : template <typename aoSystemT, class verboseT>
781 : int turbAtmosphere<aoSystemT, verboseT>::wfPS( realT ps )
782 : {
783 : static_cast<void>( ps );
784 : return 0;
785 : }
786 :
787 : template <typename aoSystemT, class verboseT>
788 : void turbAtmosphere<aoSystemT, verboseT>::nextWF( wavefront<realT> &wf )
789 : {
790 : if( !m_aosys )
791 : {
792 : throw mx::exception( error_t::paramnotset, "the AO system pointer is not set" );
793 : }
794 :
795 : if( !m_pupil )
796 : {
797 : throw mx::exception( error_t::paramnotset, "the m_pupil pointer is not set" );
798 : }
799 :
800 : shift( wf.phase, m_nWf * m_timeStep );
801 : ++m_nWf;
802 :
803 : wf.phase *= *m_pupil;
804 :
805 : realT pixVal = sqrt( m_aosys->Fg() ) * ( m_aosys->D() / m_wfSz );
806 :
807 : wf.amplitude = pixVal * ( *m_pupil );
808 : }
809 :
810 : } // namespace sim
811 : } // namespace AO
812 : } // namespace mx
813 :
814 : #endif // mx_AO_sim_turbAtmosphere_hpp
|
