Line data Source code
1 : /** \file turbLayer.hpp
2 : * \brief Declaration and definition of a turbulence layer.
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_turbLayer_hpp
30 : #define mx_AO_sim_turbLayer_hpp
31 :
32 : #include <vector>
33 :
34 : #include <Eigen/Dense>
35 :
36 : #include "../../error/error.hpp"
37 : #include "../../error/exception.hpp"
38 : #include "../../improc/imageTransforms.hpp"
39 : #include "../../math/constants.hpp"
40 : #include "../../math/func/jinc.hpp"
41 : #include "../../math/randomT.hpp"
42 : #include "../../sigproc/psdFilter.hpp"
43 : #include "../../sigproc/psdUtils.hpp"
44 :
45 : namespace mx
46 : {
47 : namespace AO
48 : {
49 : namespace sim
50 : {
51 :
52 : template <typename _aoSystemT, class _verboseT>
53 : struct turbAtmosphere;
54 :
55 : /// Simulation of a single turbulent layer
56 : /** \todo document this
57 : * \todo add facility for changing interpolator
58 : *
59 : * \ingroup mxAOSim
60 : */
61 : template <typename _aoSystemT, class _verboseT = mx::verbose::d>
62 : struct turbLayer
63 : {
64 : typedef _aoSystemT aoSystemT;
65 : typedef _verboseT verboseT;
66 :
67 : typedef typename aoSystemT::realT realT;
68 : typedef Eigen::Array<realT, -1, -1> imageT;
69 :
70 : turbAtmosphere<aoSystemT, verboseT> *m_parent{ nullptr };
71 :
72 : uint32_t m_layerNo;
73 :
74 : uint32_t m_scrnSz;
75 :
76 : imageT m_freq;
77 : imageT m_psd;
78 :
79 : realT m_dx{ 0 };
80 : realT m_dy{ 0 };
81 :
82 : realT m_x0;
83 : realT m_y0;
84 :
85 : int m_last_wdx;
86 : int m_last_wdy;
87 :
88 : imageT m_phase;
89 : imageT m_shiftPhaseWP;
90 : imageT m_shiftPhase;
91 : imageT m_shiftPhaseWork;
92 :
93 : mx::math::uniDistT<realT> uniVar; ///< Uniform deviate, used in shiftRandom.
94 :
95 : void setLayer( turbAtmosphere<aoSystemT, verboseT> *parent, int layerNo, uint32_t scrnSz );
96 :
97 : uint32_t layerNo();
98 :
99 : uint32_t scrnSz();
100 :
101 : void alloc();
102 :
103 : /// Deallocate memory necessary for phase screen generation.
104 : void genDealloc();
105 :
106 : /// Shift to a timestep.
107 : /**
108 : * \param [in] dt is the new timestep.
109 : */
110 : void shift( realT dt );
111 :
112 : /// Seed the uniform deviation. Call this if you intend to use shiftRandom.
113 : /** This only needs to be called once.
114 : */
115 : void initRandom();
116 :
117 : /// Shift by a random amount using the uniform distribution
118 : /** Call initRandom() once before calling this method.
119 : *
120 : * \param [in] nofract if true then the fractional part is ignored, and only a whole-pixel shift is executed.
121 : * Default=false
122 : */
123 : void shiftRandom( bool nofract = false );
124 : };
125 :
126 : // template<typename aoSystemT>
127 : // turbLayer<aoSystemT, verboseT>::turbLayer()
128 : // {
129 : // }
130 :
131 : template <typename aoSystemT, class verboseT>
132 2 : void turbLayer<aoSystemT, verboseT>::setLayer( turbAtmosphere<aoSystemT, verboseT> *parent,
133 : int layerNo,
134 : uint32_t scrnSz )
135 : {
136 2 : m_parent = parent;
137 2 : m_layerNo = layerNo;
138 :
139 2 : m_scrnSz = scrnSz;
140 :
141 2 : if( m_parent == nullptr )
142 : {
143 0 : throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent is nullptr)" );
144 : }
145 :
146 2 : if( m_parent->aosys() == nullptr )
147 : {
148 0 : throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent->aosys() is nullptr)" );
149 : }
150 :
151 2 : realT vwind = m_parent->aosys()->atm.layer_v_wind( m_layerNo );
152 2 : realT dwind = m_parent->aosys()->atm.layer_dir( m_layerNo );
153 2 : realT pupD = m_parent->aosys()->D();
154 2 : uint32_t wfSz = m_parent->wfSz();
155 :
156 2 : m_dx = vwind * cos( dwind ) / ( pupD / wfSz );
157 2 : m_dy = vwind * sin( dwind ) / ( pupD / wfSz );
158 :
159 2 : alloc();
160 2 : }
161 :
162 : template <typename aoSystemT, class verboseT>
163 : uint32_t turbLayer<aoSystemT, verboseT>::layerNo()
164 : {
165 : return m_layerNo;
166 : }
167 :
168 : template <typename aoSystemT, class verboseT>
169 10 : uint32_t turbLayer<aoSystemT, verboseT>::scrnSz()
170 : {
171 10 : return m_scrnSz;
172 : }
173 :
174 : template <typename aoSystemT, class verboseT>
175 2 : void turbLayer<aoSystemT, verboseT>::alloc()
176 : {
177 2 : if( m_parent == nullptr )
178 : {
179 0 : throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent is nullptr)" );
180 : }
181 :
182 2 : if( m_parent->aosys() == nullptr )
183 : {
184 0 : throw mx::exception<verboseT>( error_t::paramnotset, "parent is not set (m_parent->aosys() is nullptr)" );
185 : }
186 :
187 2 : m_phase.resize( m_scrnSz, m_scrnSz );
188 :
189 2 : uint32_t wfSz = m_parent->wfSz();
190 :
191 2 : uint32_t buffSz = m_parent->buffSz();
192 :
193 2 : m_shiftPhaseWP.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
194 :
195 2 : m_shiftPhase.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
196 :
197 2 : m_shiftPhaseWork.resize( wfSz + 2 * buffSz, wfSz + 2 * buffSz );
198 :
199 2 : initRandom();
200 :
201 2 : m_x0 = 0;
202 2 : m_y0 = 0;
203 2 : m_last_wdx = m_scrnSz + 1;
204 2 : m_last_wdy = m_scrnSz + 1;
205 :
206 : // Now set up for generation
207 2 : realT r0 = m_parent->aosys()->atm.r_0( m_parent->aosys()->lam_sci() );
208 2 : realT L0 = m_parent->aosys()->atm.L_0( m_layerNo );
209 2 : realT l0 = m_parent->aosys()->atm.l_0( m_layerNo );
210 :
211 2 : m_psd.resize( m_scrnSz, m_scrnSz );
212 :
213 2 : m_freq.resize( m_scrnSz, m_scrnSz );
214 2 : sigproc::frequencyGrid<imageT>( m_freq, m_parent->aosys()->D() / wfSz );
215 :
216 2 : realT beta = 0.0218 / pow( r0, 5. / 3. );
217 2 : realT sqrt_alpha = 0.5 * 11. / 3.;
218 :
219 : realT L02;
220 2 : if( L0 > 0 )
221 : {
222 0 : L02 = 1.0 / ( L0 * L0 );
223 : }
224 : else
225 : {
226 2 : L02 = 0;
227 : }
228 :
229 : // clang-format off
230 2 : #pragma omp parallel for // clang-format on
231 : for( size_t jj = 0; jj < m_scrnSz; ++jj )
232 : {
233 : for( size_t ii = 0; ii < m_scrnSz; ++ii )
234 : {
235 : realT p;
236 : if( m_freq( ii, jj ) == 0 && L02 == 0 )
237 : {
238 : p = 0;
239 : }
240 : else
241 : {
242 : p = beta / pow( pow( m_freq( ii, jj ), 2 ) + L02, sqrt_alpha );
243 : if( l0 > 0 )
244 : {
245 : p *= exp( -1 * pow( m_freq( ii, jj ) * l0, 2 ) );
246 : }
247 : }
248 :
249 : realT Ppiston = 0;
250 : realT Ptiptilt = 0;
251 : if( m_parent->aosys()->psd.subPiston() )
252 : {
253 : Ppiston =
254 : pow( 2 * math::func::jinc( math::pi<realT>() * m_freq( ii, jj ) * m_parent->aosys()->D() ), 2 );
255 : }
256 : if( m_parent->aosys()->psd.subTipTilt() )
257 : {
258 : Ptiptilt =
259 : pow( 4 * math::func::jincN( 2, math::pi<realT>() * m_freq( ii, jj ) * m_parent->aosys()->D() ), 2 );
260 : }
261 :
262 : m_psd( ii, jj ) = sqrt( p * ( 1.0 - Ppiston - Ptiptilt ) );
263 : }
264 : }
265 :
266 2 : if( m_parent->shLevel() > 0 )
267 : {
268 2 : realT onoff = 1;
269 2 : if( m_parent->outerSubHarmonics() )
270 : {
271 2 : onoff = 0; // just make all of these 0
272 : }
273 :
274 2 : m_psd( 0, 0 ) = 0;
275 2 : m_psd( 0, 1 ) *= 0.5 * 0.5 * onoff;
276 2 : m_psd( 1, 0 ) *= 0.5 * 0.5 * onoff;
277 2 : m_psd( m_psd.rows() - 1, 0 ) *= 0.5 * 0.5 * onoff;
278 2 : m_psd( 0, m_psd.cols() - 1 ) *= 0.5 * 0.5 * onoff;
279 :
280 2 : m_psd( 1, 1 ) *= 0.75 * 0.75 * onoff;
281 2 : m_psd( m_psd.rows() - 1, m_psd.cols() - 1 ) *= 0.75 * 0.75 * onoff;
282 2 : m_psd( m_psd.rows() - 1, 1 ) *= 0.75 * 0.75 * onoff;
283 2 : m_psd( 1, m_psd.cols() - 1 ) *= 0.75 * 0.75 * onoff;
284 : }
285 2 : }
286 :
287 : template <typename aoSystemT, class verboseT>
288 0 : void turbLayer<aoSystemT, verboseT>::genDealloc()
289 : {
290 0 : m_freq.resize( 0, 0 );
291 0 : m_psd.resize( 0, 0 );
292 0 : }
293 :
294 : template <typename aoSystemT, class verboseT>
295 : void turbLayer<aoSystemT, verboseT>::shift( realT dt )
296 : {
297 : int wdx, wdy;
298 : realT ddx, ddy;
299 :
300 : ddx = m_x0 + m_dx * dt;
301 : ddy = m_y0 + m_dy * dt;
302 :
303 : wdx = static_cast<int>( std::floor( ddx ) );
304 : ddx -= wdx;
305 :
306 : wdy = static_cast<int>( std::floor( ddy ) );
307 : ddy -= wdy;
308 :
309 : wdx %= static_cast<int>( m_scrnSz );
310 : wdy %= static_cast<int>( m_scrnSz );
311 :
312 : if( wdx < 0 )
313 : wdx += m_scrnSz;
314 : if( wdy < 0 )
315 : wdy += m_scrnSz;
316 :
317 : if( dt == 0 )
318 : {
319 : m_shiftPhase = m_phase;
320 : // Prime both the whole-pixel and shifted buffers so the first
321 : // subsequent sub-pixel shift has valid input data.
322 : improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
323 : m_shiftPhase = m_shiftPhaseWP;
324 : }
325 : else
326 : {
327 : // Check for a new whole-pixel shift
328 : if( wdx != m_last_wdx || wdy != m_last_wdy || m_last_wdx == m_scrnSz + 1 || m_last_wdy == m_scrnSz + 1 )
329 : {
330 : // Need a whole pixel shift (this also extracts the m_wfSz + m_buffSz subarray)
331 : improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
332 : }
333 :
334 : // Do the sub-pixel shift
335 : improc::imageShift( m_shiftPhase, m_shiftPhaseWP, ddx, ddy, improc::cubicConvolTransform<realT>( -0.5 ) );
336 : }
337 :
338 : m_last_wdx = wdx;
339 : m_last_wdy = wdy;
340 : }
341 :
342 : template <typename aoSystemT, class verboseT>
343 2 : void turbLayer<aoSystemT, verboseT>::initRandom()
344 : {
345 2 : uniVar.seed();
346 2 : }
347 :
348 : template <typename aoSystemT, class verboseT>
349 : void turbLayer<aoSystemT, verboseT>::shiftRandom( bool nofract )
350 : {
351 : int wdx, wdy;
352 : realT ddx, ddy;
353 :
354 : ddx = uniVar * ( m_scrnSz );
355 : ddy = uniVar * ( m_scrnSz );
356 :
357 : wdx = (int)trunc( ddx );
358 : ddx -= wdx;
359 : wdy = (int)trunc( ddy );
360 : ddy -= wdy;
361 :
362 : if( nofract )
363 : {
364 : ddx = 0;
365 : ddy = 0;
366 : }
367 :
368 : improc::imageShiftWP( m_shiftPhaseWP, m_phase, wdx, wdy );
369 :
370 : if( ddx != 0 || ddy != 0 )
371 : {
372 : improc::imageShift( m_shiftPhase, m_shiftPhaseWP, ddx, ddy, improc::cubicConvolTransform<realT>( -0.5 ) );
373 : }
374 : else
375 : {
376 : m_shiftPhase = m_shiftPhaseWP;
377 : }
378 : }
379 :
380 : } // namespace sim
381 : } // namespace AO
382 : } // namespace mx
383 :
384 : #endif // mx_AO_sim_turbLayer_hpp
|
