Line data Source code
1 : /** \file aoSystem.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Declares and defines an analytical AO system
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015-2022 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef aoSystem_hpp
28 : #define aoSystem_hpp
29 :
30 : #include "../../mxlib.hpp"
31 : #include "../../math/constants.hpp"
32 : #include "../../math/roots.hpp"
33 :
34 : #include "aoConstants.hpp"
35 :
36 : #include "aoAtmosphere.hpp"
37 : #include "aoPSDs.hpp"
38 : #include "aoWFS.hpp"
39 :
40 : #define FITTING_ERROR_NO 0
41 : #define FITTING_ERROR_ZERO 1
42 : #define FITTING_ERROR_X 2
43 : #define FITTING_ERROR_Y 3
44 :
45 : namespace mx
46 : {
47 : namespace AO
48 : {
49 : namespace analysis
50 : {
51 :
52 : /// Describes an analytic adaptive optics (AO) system.
53 : /**
54 : * Templatized by the turbulence power spectral density (PSD).
55 : *
56 : * \tparam realT the floating point type used for all calculations
57 : * \tparam inputSpecT specifies the turbulence spatial PSD type
58 : * \tparam iosT is an output stream type with operator \<\< defined (default is std::ostream)
59 : *
60 : * \ingroup mxAOAnalytic
61 : */
62 : template <typename _realT, class _inputSpectT, typename iosT = std::ostream>
63 : class aoSystem
64 : {
65 : typedef verbose::d verboseT;
66 :
67 : public:
68 : typedef _realT realT;
69 : typedef _inputSpectT inputSpectT;
70 : typedef aoAtmosphere<realT> aoAtmosphereT;
71 :
72 : aoAtmosphereT atm;
73 : inputSpectT psd;
74 :
75 : protected:
76 : realT m_D{ 0 }; ///< Telescope diameter [m]
77 :
78 : std::vector<realT> m_d_min{ 0 }; ///< Minimum AO system actuator pitch [m]. One per WFS mode.
79 :
80 : realT m_d_opt{ 0 }; ///< Current optimum AO system actuator pitch [m]
81 : bool m_optd{ false }; ///< Flag controlling whether actuator pitch is optimized (true) or just uses m_d_min (false).
82 : ///< Default: true.
83 : realT m_optd_delta{ 1.0 }; ///< The fractional change from d_min used in optimization. Set to 1 for integer
84 : ///< binnings, > 1 for finer sampling.
85 :
86 : wfs<realT, iosT> *m_wfsBeta{ nullptr }; ///< The WFS beta_p class.
87 : bool m_ownWfsBeta{ false }; ///< Flag indicating if the WFS beta_p pointer is owned by this instance.
88 :
89 : realT m_opticalGain{ 1 }; /**< The optical gain of the WFS, 0-1. Treated as a sensitivity
90 : reduction in measurement noise. Default 1.*/
91 :
92 : protected:
93 : realT m_lam_wfs{ 0 }; ///< WFS wavelength [m]
94 :
95 : // WFS detector parameters, which may be a function of binning
96 : std::vector<realT> m_npix_wfs{ 0 }; ///< Number of WFS pixels. One per WFS mode.
97 : std::vector<realT> m_ron_wfs{ 0 }; ///< WFS readout noise [electrons/pix]. One per WFS mode.
98 : std::vector<realT> m_Fbg{ 0 }; ///< Background flux, [photons/sec/pixel]. One per WFS mode.
99 : std::vector<realT> m_minTauWFS{ 0 }; ///< Minimum WFS exposure time [sec]. One per WFS mode.
100 :
101 : bool m_bin_npix{ true }; ///< Flag controlling whether or not to bin WFS pixels according to the actuator spacing.
102 :
103 : int m_bin_opt{ 0 }; /**< The optimum binning factor. If WFS modes are used, this is the mode
104 : index (0 to N-1). If not, it is 1 minus the pixel binning factor.
105 : It is always 1 minus the actuator binning factor. */
106 :
107 : realT m_tauWFS{ 0 }; ///< Actual WFS exposure time [sec]
108 :
109 : realT m_deltaTau{ 0 }; ///< Loop latency [sec]
110 :
111 : bool m_optTau{ true }; ///< Flag controlling whether optimum integration time is calculated (true) enforcing
112 : ///< m_minTauWFS, or if m_tauWFS is used (false). Default: true.
113 :
114 : realT m_lam_sci{ 0 }; ///< Science wavelength [m]
115 :
116 : realT m_zeta{ 0 }; ///< Zenith angle [radians]
117 : realT m_secZeta{ 1 }; ///< Secant of the Zenith angle (calculated)
118 :
119 : int m_fit_mn_max{ 100 }; ///< Maximum spatial frequency index to use for fitting error calculation.
120 :
121 : bool m_circularLimit{ false }; ///< Flag to indicate that the spatial frequency limit is circular, not square.
122 :
123 : realT m_spatialFilter_ku{ std::numeric_limits<realT>::max() }; ///< The spatial filter cutoff in u, [m^-1]
124 : realT m_spatialFilter_kv{ std::numeric_limits<realT>::max() }; ///< The spatial filter cutoff in v, [m^-1]
125 :
126 : realT m_ncp_wfe{ 0 }; ///< Static WFE [m rms]
127 : realT m_ncp_alpha{ 2 }; ///< Power-law exponent for the NCP aberations. Default is 2.
128 :
129 : realT m_F0{ 0 }; ///< 0 mag flux from star at WFS [photons/sec]
130 :
131 : realT m_starMag{ 0 }; ///< The magnitude of the star.
132 :
133 : bool m_specsChanged{ true }; ///< Flag to indicate that a specification has changed.
134 : bool m_dminChanged{ true }; ///< Flag to indicate that d_min has changed.
135 :
136 : Eigen::Array<int, -1, -1> m_controlledModes; ///< Map of which modes are under control. Set by calcStrehl.
137 :
138 : realT m_wfeMeasurement{ 0 }; ///< Total WFE due to measurement a error [rad^2 at m_lam_sci]
139 : realT m_wfeTimeDelay{ 0 }; ///< Total WFE due to time delay [rad^2 at m_lam_sci]
140 : realT m_wfeFitting{ 0 }; ///< Total WFE due to fitting error [rad^2 at m_lam_sci]
141 : realT m_wfeChromScintOPD{ 0 }; ///< Total WFE due to the chromaticity of scintillation OPD [rad^2 at lam_sc]
142 : realT m_wfeChromIndex{ 0 }; ///< Total WFE due to the chromaticity of the index of refraction [rad^2 at lam_Sci]
143 : realT m_wfeAnisoOPD{ 0 }; ///< Total WFE due to dispersive anisoplanatism OPD.
144 :
145 : realT m_wfeNCP{ 0 }; ///< Total WFE due to NCP errors [rad^2 at m_lam_sci]
146 :
147 : realT m_wfeVar{ 0 }; ///< The WFE variance, in meters^2. Never use this directly, instead use wfeVar().
148 : realT m_strehl{ 0 }; ///< Strehl ratio, a calculated quantity. Never use this directdy, instead use strehl().
149 :
150 : public:
151 : /// Default c'tor
152 : /** Calls initialize().
153 : */
154 : aoSystem();
155 :
156 : /// Destructor
157 : ~aoSystem();
158 :
159 : /// Load the default parameters from Guyon, 2005 \cite guyon_2005.
160 : /**
161 : *
162 : */
163 : void loadGuyon2005();
164 :
165 : /// Load parameters corresponding to the MagAO-X system.
166 : void loadMagAOX();
167 :
168 : /// Load parameters corresponding to the G-MagAO-X system.
169 : void loadGMagAOX();
170 :
171 : /// Set the value of the primary mirror diameter.
172 : /**
173 : */
174 : void D( realT nD /**< [in] is the new value of m_D. */ );
175 :
176 : /// Get the value of the primary mirror diamter
177 : /**
178 : * \returns the current value of m_D.
179 : */
180 : realT D();
181 :
182 : /// Set the minimum subaperture sampling for each WFS mode.
183 : /** This sets the vector
184 : */
185 : void d_min( const std::vector<realT> &nd /**< [in] is the new values of m_d_min */ );
186 :
187 : /// Set the minimum subaperture sampling for a given WFS mode
188 : /** This sets the entry in the vector, if it has the proper length.
189 : */
190 : void d_min( int idx, ///< [in] the index of the WFS mode
191 : realT nd ///< [in] is the new value of m_d_min[idx]
192 : );
193 :
194 : /// Get the minimum subaperture sampling for a given WFS mode
195 : /**
196 : * \returns the current value of m_d_min[idx]
197 : */
198 : realT d_min( size_t idx /**< [in] the index of the WFS mode.*/ );
199 :
200 : /// Get the minimum subaperture sampling for all WFS modes.
201 : /**
202 : * \returns the current value of m_d_min.
203 : */
204 : std::vector<realT> d_min();
205 :
206 : /// Set whether or not the value of d is optimized or just set to m_d_min.
207 : /**
208 : */
209 : void optd( bool od /**< [in] is the new value of m_optd */ );
210 :
211 : /// Get the value of m_optd.
212 : /**
213 : * \returns the new value of m_optd_delta.
214 : */
215 : bool optd();
216 :
217 : /// Set the fractional change in actuator spacing for optimization.
218 : /** Sets the fraction of m_d_min by which the optimizer changes actautor spacing.
219 : */
220 : void optd_delta( realT odd /**< [in] is the new value of m_optd_delta */ );
221 :
222 : /// Get the value of the fractional change in actuator spacing for optimization..
223 : /**
224 : * \returns the value of m_optd_delta.
225 : */
226 : realT optd_delta();
227 :
228 : /// Set the WFS Beta pointer
229 : /** The WFS beta parameter determines the photon noise sensitivity.
230 : *
231 : * \tparam wfsT is a type derived from ao::analysis::wfs
232 : */
233 : template <typename wfsT>
234 : void wfsBeta( const wfsT &w /**< [in] an object derived from ao::analysis::wfs base class*/ );
235 :
236 : /// Set the WFS Beta pointer
237 : /** The WFS beta parameter determines the photon noise sensitivity.
238 : *
239 : * \tparam wfsT is a type derived from ao::analysis::wfs
240 : */
241 : template<typename wfsT>
242 : void wfsBeta( const wfsT * w /**< [in] pointer to an object derived from ao::analysis::wfs base class. If nullptr then a new object is allocated and managed.*/);
243 :
244 : /// Get the WFS Beta pointer
245 : /**
246 : * \returns a pointer to the current m_wfsBeta
247 : */
248 : wfs<realT, iosT> *wfsBeta();
249 :
250 : /// Get the value of beta_p for a spatial frequency
251 : /** beta_p is the photon noise sensitivity of the WFS.
252 : *
253 : * \returns beta_p as calculated by the WFS.
254 : */
255 : realT beta_p( realT m, ///< [in] the spatial frequency index
256 : realT n ///< [in] the spatial frequency index
257 : );
258 :
259 : /// Get the value of beta_r for a spatial frequency
260 : /** beta_r is the read noise sensitivity of the WFS.
261 : *
262 : * \returns beta_r as calculated by the WFS.
263 : */
264 : realT beta_r( realT m, ///< [in] the spatial frequency index
265 : realT n ///< [in] the spatial frequency index
266 : );
267 :
268 : /// Set the value of the Optical Gain.
269 : /**
270 : */
271 : void opticalGain( realT og /**< [in] the new value of m_opticalGain*/ );
272 :
273 : /// Get the value of the optical gain
274 : /** Optical gain is applied in the WFS noise calculation
275 : *
276 : * \returns the current value of m_opticalGain
277 : */
278 : realT opticalGain();
279 :
280 : /// Set the value of the WFS wavelength.
281 : /**
282 : */
283 : void lam_wfs( realT nlam /**< [in] is the new value of m_lam_wfs */ );
284 :
285 : /// Get the value of the WFS wavelength.
286 : /**
287 : * \returns the current value of m_lam_wfs.
288 : */
289 : realT lam_wfs();
290 :
291 : /// Set the number of pixels in the WFS for each WFS mode
292 : /** This sets the vector.
293 : */
294 : void npix_wfs( const std::vector<realT> &npix /**< [in] is the new values of m_npix_wfs */ );
295 :
296 : /// Set the number of pixels in the WFS for a given WFS mode
297 : /** This sets the entry in the vector, if it has the proper length.
298 : */
299 : void npix_wfs( int idx, ///< [in] the index of the WFS mode
300 : realT npix ///< [in] is the new value of m_npix_wfs[idx]
301 : );
302 :
303 : /// Get the number of pixels in the WFS for a given WFS mode
304 : /**
305 : * \returns the current value of m_npix_wfs[idx]
306 : */
307 : realT npix_wfs( size_t idx /**< [in] the index of the WFS mode.*/ );
308 :
309 : /// Get the number of pixels in the WFS for all WFS modes
310 : /**
311 : * \returns the current value of m_npix_wfs
312 : */
313 : std::vector<realT> npix_wfs();
314 :
315 : /// Set the value of the WFS readout noise for each WFS mode
316 : /** This sets the vector.
317 : */
318 : void ron_wfs( const std::vector<realT> &nron /**< [in] is the new value of m_ron_wfs */ );
319 :
320 : /// Set the value of the WFS readout noise for a given bining mode
321 : /** This sets the entry in the vector if the vector has the proper length
322 : */
323 : void ron_wfs( int idx, ///< [in] the index of the WFS mode
324 : realT nron ///< [in] is the new value of m_ron_wfs [idx]
325 : );
326 :
327 : /// Get the value of the WFS readout noise for a given WFS mode
328 : /**
329 : * \returns the current value of m_ron_wfs[idx]
330 : */
331 : realT ron_wfs( size_t idx /**< [in] the index of the WFS mode.*/ );
332 :
333 : /// Get the value of the WFS readout noise for all WFS modes
334 : /**
335 : * \returns the current value of m_ron_wfs
336 : */
337 : std::vector<realT> ron_wfs();
338 :
339 : /// Set the value of the background fluxes.
340 : /** This sets the value of the background flux for each WFS mode
341 : */
342 : void Fbg( const std::vector<realT> &fbg /**< [in] is the new value of m_Fbg */ );
343 :
344 : /// Set a single value of the background flux for a given WFS mode
345 : /** This sets the entry in the vector if the vector has the proper length
346 : */
347 : void Fbg( int idx, ///< [in] the index of the WFS mode
348 : realT fbg ///< [in] is the new value of m_Fbg[idx]
349 : );
350 :
351 : /// Get the value of the background flux for a given WFS mode
352 : /**
353 : * \returns m_Fbg[idx]
354 : */
355 : realT Fbg( size_t idx /**< [in] the index of the WFS mode.*/ );
356 :
357 : /// Get the value of the background flux for all WFS modes
358 : /**
359 : * \returns the current value of m_Fbg
360 : */
361 : std::vector<realT> Fbg();
362 :
363 : /// Set the value of the minimum WFS exposure times.
364 : /** This sets the vector of the minimum WFS exposure times
365 : */
366 : void minTauWFS( const std::vector<realT> &ntau /**< [in] is the new value of m_minTauWFS */ );
367 :
368 : /// Set a single value of the minimum WFS exposure time for a given WFS mode.
369 : /** This sets the entry in the vector if the vector has sufficient length.
370 : */
371 : void minTauWFS( size_t idx, ///< [in] the index of the WFS mode
372 : realT ntau ///< [in] is the new value of m_minTauWFS
373 : );
374 :
375 : /// Get the value of the minimum WFS exposure time for a given binning.
376 : /**
377 : * \returns the current value of m_minTauWFS[idx].
378 : */
379 : realT minTauWFS( size_t idx /**< [in] the index of the WFS mode.*/ );
380 :
381 : /// Get the values of the minimum WFS exposure time.
382 : /**
383 : * \returns the current value of m_minTauWFS.
384 : */
385 : std::vector<realT> minTauWFS();
386 :
387 : /// Set the value of the pixel binning flag
388 : /**
389 : */
390 : void bin_npix( bool bnp /**< [in] is the new value of m_bin_npix */ );
391 :
392 : /// Get the value of the pixel binngin flag
393 : /**
394 : * \returns
395 : */
396 : bool bin_npix();
397 :
398 : /// Get the value of the optimum binning factor/index.
399 : /** If WFS modes are used, this is the mode index (0 to N-1). If not, it is 1 minus the pixel binning factor.
400 : * It is always 1 minus the actuator binning factor.
401 : * This calls opt_d() to perform the optimization if needed.
402 : *
403 : * \returns the current value of m_bin_opt.
404 : */
405 : int bin_opt();
406 :
407 : /// Set the value of the WFS exposure time.
408 : /**
409 : */
410 : void tauWFS( realT ntau /**< [in] is the new value of m_tauWFS */ );
411 :
412 : /// Get the value of the minimum WFS exposure time.
413 : /**
414 : * \returns the current value of m_tauWFS.
415 : */
416 : realT tauWFS();
417 :
418 : /// Set the value of m_deltaTau.
419 : /**
420 : */
421 : void deltaTau( realT ndel /**< [in] is the new value of m_deltaTau*/ );
422 :
423 : /// Get the value of m_deltaTau.
424 : /**
425 : * \returns the current value of m_deltaTau.
426 : */
427 : realT deltaTau();
428 :
429 : /// Set the value of m_optTau.
430 : /**
431 : */
432 : void optTau( bool ot /**< [in] is the new value of m_optTau */ );
433 :
434 : /// Get the value of m_optTau.
435 : /**
436 : * \returns the current value of m_optTau.
437 : */
438 : bool optTau();
439 :
440 : /// Set the science wavelength.
441 : /**
442 : */
443 : void lam_sci( realT nlam /**< [in] is the new value of m_lam_sci */ );
444 :
445 : /// Get the science wavelength.
446 : /**
447 : * \returns the current value of m_lam_sci.
448 : */
449 : realT lam_sci();
450 :
451 : /// Set the zenith angle, and its secant.
452 : /**
453 : */
454 : void zeta( realT nz /**< [in] The new value of m_zeta */ );
455 :
456 : /// Get the zenith angle
457 : /**
458 : * \return the current value of m_zeta
459 : */
460 : realT zeta();
461 :
462 : /// Get the zecant of the zenith angle
463 : /**
464 : * \return the current value of m_secZeta
465 : */
466 : realT secZeta();
467 :
468 : /// Set the value of m_fit_mn_max
469 : /**
470 : */
471 : void fit_mn_max( int mnm /**< [in] is the new value of m_fit_mn_max */ );
472 :
473 : /// Get the value of m_fit_mn_max
474 : /**
475 : */
476 : int fit_mn_max();
477 :
478 : /// Set the value of the circularLimit flag
479 : /** If this is true, then the spatial frequency limits are treated as a circle
480 : * rather than a square. This will create a circular dark hole.
481 : */
482 : void circularLimit( bool cl /**< [in] the new value of the circularLimit flag*/ );
483 :
484 : /// Get the value of the circularLimit flag
485 : /**
486 : * \returns the current value of m_circularLimit
487 : */
488 : bool circularLimit();
489 :
490 : /// Set the value of spatialFilter_ku
491 : /**
492 : */
493 : void spatialFilter_ku( realT ku /**< [in] is the new value of spatialFilter_ku*/ );
494 :
495 : /// Get the value of spatialFilter_ku
496 : /** \returns the current value of m_spatialFilter_ku
497 : */
498 : realT spatialFilter_ku();
499 :
500 : /// Set the value of spatialFilter_kv
501 : /**
502 : */
503 : void spatialFilter_kv( realT kv /**< [in] is the new value of spatialFilter_kv*/ );
504 :
505 : /// Get the value of spatialFilter_kv
506 : /** \returns the current value of m_spatialFilter_kv
507 : */
508 : realT spatialFilter_kv();
509 :
510 : /// Set the value of the non-common path WFE.
511 : /**
512 : */
513 : void ncp_wfe( realT nwfe /**< [in] is the new value of m_ncp_wfe*/ );
514 :
515 : /// Get the value of the non-common path WFE.
516 : /**
517 : * \returns the current value of m_ncp_wfe.
518 : */
519 : realT ncp_wfe();
520 :
521 : /// Set the value of the non-common path WFE PSD index.
522 : /**
523 : */
524 : void ncp_alpha( realT alpha /**< [in] is the new value of m_ncp_alpha*/ );
525 :
526 : /// Get the value of the non-common path WFE PSD index.
527 : /**
528 : * \returns the current value of m_ncp_alpha.
529 : */
530 : realT ncp_alpha();
531 :
532 : /// Set the value of the 0 magnitude photon rate
533 : /** This is the photon rate at the WFS, photons/sec.
534 : *
535 : */
536 : void F0( realT nF0 /**< [in] is the new value of m_F0.*/ );
537 :
538 : /// Get the value of the 0 magnitude photon rate
539 : /**
540 : * \returns the current value of m_F0.
541 : */
542 : realT F0();
543 :
544 : /// Set the value of the Star's magnitude
545 : /**
546 : */
547 : void starMag( realT nmag /**< [in] is the new value of m_starMag.*/ );
548 :
549 : /// Get the value of the Star's magnitude
550 : /**
551 : * \returns the current value of m_starMag
552 : */
553 : realT starMag();
554 :
555 : /// The photon flux at a given star magnitude.
556 : /**
557 : * \returns the photon flux of the star in the bandpass assumed by F0.
558 : */
559 : realT Fg( realT mag /**< [in] is the magnitude of the star. */ );
560 :
561 : /// Get the photon rate at the current Star magnitude.
562 : /** Calculates \f$ F_\gamma = F_0 10^{-0.4 m} \f$ where \f$ F_0 \f$ is m_F0 and \f$ m \f$ is m_starMag.
563 : *
564 : * \returns the current value of the current photon rate.
565 : */
566 : realT Fg();
567 :
568 : /// Access the array of controlledModes
569 : /** This array indicates which modes are controlled after optimizeStrehl chooses them.
570 : *
571 : * \returns a const reference to m_controlledModes.
572 : */
573 : const Eigen::Array<int, -1, -1> &controlledModes();
574 :
575 : /** \name Measurement Error
576 : * Calculating the WFE due to WFS measurement noise.
577 : * @{
578 : */
579 :
580 : /// Calculate the terms of the signal to noise ratio squared (S/N)^2 for the WFS measurement
581 : /** The S/N squared is
582 : \f[
583 : (S/N)^2 = \frac{ F_\gamma^2 \tau_{wfs}^2 }{ F_\gamma \tau_{wfs} + n_{pix} F_{bg} \tau_{wfs} + n_{pix}
584 : \sigma_{ron}^2 } \f]
585 :
586 : *
587 : * \returns the S/N squared
588 : */
589 : realT
590 : signal2Noise2( realT &Nph, ///< [out] the number of photons
591 : realT tau_wfs, ///< [in] specifies the WFS exposure time. If 0, then optimumTauWFS is used
592 : realT d, ///< [in] the actuator spacing in meters, used if binning WFS pixels
593 : int b ///< [in] the binning parameter. Either the WFS mode index, or the binning factor minus 1.
594 : );
595 :
596 : /// Calculate the measurement noise at a spatial frequency and specified actuator spacing
597 : /** Calculates the wavefront phase variance due measurement noise at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
598 : *
599 : * \returns the measurement noise in rad^2 rms at the science wavelength
600 : */
601 : realT
602 : measurementError( realT m, ///< [in] specifies the u component of the spatial frequency
603 : realT n, ///< [in] specifies the v component of the spatial frequency
604 : realT d, ///< [in] the actuator spacing in meters
605 : int b ///< [in] the binning parameter. Either the WFS mode index, or the binning factor minus 1.
606 : );
607 :
608 : /// Calculate the total measurement error over all corrected spatial frequencies
609 : /** This totals the measurement WFE at the optimum actuator spacing and binning.
610 : *
611 : * \overload
612 : *
613 : * \returns the total WFE due to measurement error.
614 : */
615 : realT measurementErrorTotal();
616 :
617 : ///@}
618 :
619 : /** \name Time Delay Error
620 : * Calculating the WFE due to time delay.
621 : * @{
622 : */
623 :
624 : /// Calculate the time delay at a spatial frequency at the optimum exposure time and the specified actuator spacing
625 : /** Calculates the wavefront phase variance due to time delay at \f$ f = (m/D)\hat{u} + (n/D)\hat{v} \f$.
626 : *
627 : * \returns the measurement noise in rad^2 rms at the science wavelength
628 : */
629 : realT
630 : timeDelayError( realT m, ///< [in] specifies the u component of the spatial frequency
631 : realT n, ///< [in] specifies the v component of the spatial frequency
632 : realT d, ///< [in] the actuator spacing, in meters
633 : int b ///< [in] the binning parameter. Either the WFS mode index, or the binning factor minus 1.
634 : );
635 :
636 : /// Calculate the time delay error over all corrected spatial frequencies
637 : /** This totals the time delay WFE at the optimum actuator spacing and binning.
638 : *
639 : * \overload
640 : *
641 : * \returns the total WFE due to time delay.
642 : */
643 : realT timeDelayErrorTotal();
644 :
645 : ///@}
646 :
647 : /** \name Fitting Error
648 : * Calculating the WFE due to uncorrected spatial frequencies.
649 : * @{
650 : */
651 :
652 : /// Calculate the fitting error at a specific spatial frequency.
653 : /**
654 : * \returns the fitting error in rad^2 rms at the science wavelength at (m,n).
655 : */
656 : realT fittingError( realT m, ///< [in] specifies the u component of the spatial frequency
657 : realT n ///< [in] specifies the v component of the spatial frequency
658 : );
659 :
660 : /// Calculate the total fitting error over all uncorrected spatial frequencies.
661 : /** This totals the fitting WFE for the optimum actuator spacing and binning.
662 : *
663 : * \overload
664 : *
665 : * \returns the total fitting error.
666 : */
667 : realT fittingErrorTotal();
668 :
669 : ///@}
670 :
671 : /** \name Chromatic Errors
672 : * Calculating the WFE due to chromaticity in scintillaion, index of refraction, and dispersion.
673 : * @{
674 : */
675 :
676 : /// Calculate the wavefront error due to scintillation chromaticity in the OPD at a spatial frequency.
677 : /** This totals the WFE from scintillation chromaticity in the OPD.
678 : *
679 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
680 : */
681 : realT chromScintOPDError( int m, ///< [in] specifies the u component of the spatial frequency
682 : int n ///< [in] specifies the v component of the spatial frequency
683 : );
684 :
685 : /// Calculate the wavefront error due to scintillation chromaticity in the OPD over all spatial frequencies.
686 : /** This totals the WFE from scintillation chromaticity in the OPD for the optimum actuator spacing.
687 : *
688 : * \overload
689 : *
690 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
691 : */
692 : realT chromScintOPDErrorTotal();
693 :
694 : /// Calculate the wavefront error due to scintillation chromaticity in amplitude for a given spatial frequency.
695 : /**
696 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
697 : */
698 : realT chromScintAmpError( int m, ///< [in] specifies the u component of the spatial frequency
699 : int n ///< [in] specifies the v component of the spatial frequency
700 : );
701 :
702 : /// Calculate the wavefront error due to scintillation chromaticity in amplitude over all spatial frequencies.
703 : /**
704 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
705 : */
706 : realT chromScintAmpError();
707 :
708 : /// Calculate the wavefront error due to chromaticity in the index of refraction at a given spatial frequency.
709 : /** This totals the WFE from chromaticity in the index of refraction.
710 : *
711 : * \overload
712 : *
713 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
714 : */
715 : realT chromIndexError( int m, ///< [in] specifies the u component of the spatial frequency
716 : int n ///< [in] specifies the v component of the spatial frequency
717 : );
718 :
719 : /// Calculate the wavefront error due to chromaticity in the index of refraction at a specific spatial frequency.
720 : /** This totals the WFE from chromaticity in the index of refraction for the optimum actuator spacing.
721 : *
722 : * \overload
723 : *
724 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
725 : */
726 : realT chromIndexErrorTotal();
727 :
728 : /// Calculate the wavefront error due to dispersive anisoplanatism in the OPD at a given spatial frequency
729 : /** This calculates the WFE from dispersive anisoplanatism in the OPD.
730 : *
731 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
732 : */
733 : realT dispAnisoOPDError( int m, ///< [in] specifies the u component of the spatial frequency
734 : int n ///< [in] specifies the v component of the spatial frequency
735 : );
736 :
737 : /// Calculate the wavefront error due to dispersive anisoplanatism in the OPD over all specific spatial frequencies.
738 : /** This totals the WFE from dispersive anisoplanatism in the OPD for the optimum actuator spacing.
739 : *
740 : * \overload
741 : *
742 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
743 : */
744 : realT dispAnisoOPDErrorTotal();
745 :
746 : /// Calculate the wavefront error due to dispersive anisoplanatism in the amplitude at a specific spatial frequency.
747 : /**
748 : * \returns the WFE in rad^2 rms at the science wavelength at (m,n).
749 : */
750 : realT dispAnisoAmpError();
751 :
752 : ///@}
753 :
754 : /** \name Optimum Parameters
755 : * Functions to calculate the optimum integration time and actuator spacing.
756 : *
757 : * @{
758 : */
759 :
760 : /// Calculate the optimum exposure time for a given spatial frequency at a specified actuator spacing
761 : /** Finds the optimum exposure time at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
762 : *
763 : * \todo Check inclusion of X in parameters
764 : *
765 : * \returns the optimum expsoure time.
766 : */
767 : realT
768 : optimumTauWFS( realT m, ///< [in] is the spatial frequency index in u
769 : realT n, ///< [in] is the spatial frequency index in v
770 : realT d, ///< [in] the actuator spacing, in meters
771 : int b ///< [in] the binning parameter. Either the WFS mode index, or the binning factor minus 1.
772 : );
773 :
774 : /// Calculate the optimum exposure time for a given spatial frequency at the optimum actuator spacing.
775 : /** Finds the optimum exposure time at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
776 : * If optd == false this uses the minimum actuator spacing.
777 : *
778 : * \todo Check inclusion of X in parameters
779 : *
780 : * \returns the optimum expsoure time.
781 : */
782 : realT optimumTauWFS( realT m, ///< [in] is the spatial frequency index in u
783 : realT n ///< [in] is the spatial frequency index in v
784 : );
785 :
786 : /// Calculate the optimum actuator spacing.
787 : /** Finds the value of m_d_opt where the fitting error is less than than the combined time delay and measurement
788 : * error.
789 : *
790 : * \returns the current value of m_d_opt.
791 : */
792 : realT d_opt();
793 :
794 : /// @}
795 :
796 : /// Calculate the NCP variance at a spatial frequency.
797 : /** Finds the NCP variance at \f$ k = (m/D)\hat{u} + (n/D)\hat{v} \f$.
798 : *
799 : * \returns the NCP variance at k.
800 : */
801 : realT ncpError( int m, ///< [in] is the spatial frequency index in u
802 : int n ///< [in] is the spatial frequency index in v
803 : );
804 :
805 : /// Calculate the total NCP variance in rad^2.
806 : /**
807 : * \returns the total NCP variance.
808 : */
809 : realT ncpError();
810 :
811 : /** \name Overall WFE and Strehl Ratio
812 : * Functions to calculate the total WFE and Strehl ratio.
813 : *
814 : * @{
815 : */
816 :
817 : protected:
818 : /// Calculate the component WFE, total WFE, and Strehl ratio.
819 : /** This should only be called when something changes.
820 : */
821 : void calcStrehl();
822 :
823 : public:
824 : /// Get the current value of the total WFE variance.
825 : /** If no changes, merely returns m_wfeVar. Calls calcStrehl if there are changes.
826 : *
827 : * \returns the current value of m_wfeVar;
828 : */
829 : realT wfeVar();
830 :
831 : /// Get the Strehl ratio for a given actuator pitch and WFS WFS mode.
832 : /**
833 : * Strehl is calculated using the extended Marechal approximation:
834 : *
835 : \f[
836 : S = e^{-\sigma_{wfe}^2}
837 : \f]
838 : * where \f$ \sigma_{wfe}^2 \f$ is the WFE for the spacing and binning
839 : *
840 : * \returns the strehl for these values
841 : */
842 : realT strehl( realT d, ///< [in] the actuator spacing, in meters
843 : int b ///< [in] the binning parameter. Either the WFS mode index, or the binning factor minus 1.
844 : );
845 :
846 : /// Get the current Strehl ratio.
847 : /** If no changes, merely returns m_strehl. Calls calcStrehl if there are changes.
848 : * Strehl is calculated using the extended Marechal approximation:
849 : *
850 : \f[
851 : S = e^{-\sigma_{wfe}^2}
852 : \f]
853 : * where \f$ \sigma_{wfe}^2 \f$ is the current value of m_wfeVar.
854 : *
855 : * \returns the current value of m_strehl.
856 : */
857 : realT strehl();
858 :
859 : /// @}
860 :
861 : /** \name Contrast
862 : * Calculating contrast
863 : *
864 : * @{
865 : */
866 :
867 : /// Worker function for raw contrast fuctions
868 : /** The logic in this calculation is the same regardless of component, except for the calculation of variance.
869 : * The varFunc function pointer is used to calculate the variance, otherwise This handles the other details such
870 : * as Strehl normalization, fitting error (if appropriate), and bounds checking.
871 : *
872 : * \note this gives the raw PSD, it must be convolved with the PSF for the true contrast.
873 : *
874 : * \tparam varFuncT is a function-pointer-to-member (of this class) with signature realT(*f)(realT, realT)
875 : */
876 : template <typename varFuncT>
877 : realT C_( int m, ///< [in] is the spatial frequency index in u/
878 : int n, ///< [in] is the spatial frequency index in v.
879 : bool normStrehl, ///< [in] flag controls whether the contrast is normalized by Strehl ratio/
880 : varFuncT varFunc, ///< [in] the variance function to use.
881 : int doFittingError ///< [in] flag to describe how fitting error should be considered for this term:
882 : ///< FITTING_ERROR_NO, FITTING_ERROR_ZERO, FITTING_ERROR_X, or FITTING_ERROR_Y.
883 : );
884 :
885 : /// Worker function for the contrast-map functions.
886 : /** The map calculation is the same for all terms, except for the calculation of variance at each pixel.
887 : * The Cfunc function pointer is used to actually get the contrast.
888 : *
889 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
890 : *
891 : * \tparam imageT is an Eigen-like image
892 : * \tparam CfuncT is a function-pointer-to-member (of this class) with signature realT (*f)(realT, realT, bool)
893 : */
894 : template <typename imageT, typename CfuncT>
895 : void C_Map( imageT &im, ///< [out] the map image to be filled in.
896 : CfuncT Cfunc, ///< [in] the raw contrast function to use for filling in the map.
897 : bool normStrehl ///< [in] flag controlling whether or not the map is normalized by Strehl
898 : );
899 :
900 : /// Calculate the residual variance due to uncorrected phase at a spatial frequency.
901 : /** Used to calculate contrast \ref C0().
902 : *
903 : * \returns variance at (m,n).
904 : */
905 : realT C0var( realT m, ///< [in] is the spatial frequency index in u
906 : realT n ///< [in] is the spatial frequency index in v
907 : );
908 :
909 : /// Calculate the contrast due to uncorrected phase, C0.
910 : /** Contrast C0 is the uncorrected phase, with the effects of scintillation included. See Guyon (2005)
911 : * \cite guyon_2005 and the updated derivation in Males \& Guyon (2018) \cite males_guyon_2018
912 : *
913 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
914 : *
915 : * \returns C0.
916 : */
917 : realT C0( realT m, ///< [in] is the spatial frequency index in u
918 : realT n, ///< [in] is the spatial frequency index in v
919 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio
920 : );
921 :
922 : /// Calculate a 2D map of contrast C0
923 : /**
924 : * \tparam imageT is an Eigen-like image type.
925 : */
926 : template <typename imageT>
927 : void C0Map( imageT &map, ///< [in] the map image to be filled in with contrast
928 : bool normStrehl );
929 :
930 : /// Calculate the residual variance due to uncorrected amplitude at a spatial frequency.
931 : /** Used to calculate contrast \ref C1().
932 : *
933 : * \returns variance at (m,n).
934 : */
935 : realT C1var( realT m, ///< [in] is the spatial frequency index in u
936 : realT n ///< [in] is the spatial frequency index in v
937 : );
938 :
939 : /// Calculate the contrast due to uncorrected amplitude, C1.
940 : /** Contrast C1 is the uncorrected amplitude due to scintillation. See Guyon (2005) \cite guyon_2005, and the
941 : * updated derivation in Males \& Guyon (2018) \cite males_guyon_2018.
942 : *
943 : * \returns C0.
944 : */
945 : realT C1( realT m, ///< [in] is the spatial frequency index in u
946 : realT n, ///< [in] is the spatial frequency index in v
947 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio
948 : );
949 :
950 : /// Calculate a 2D map of contrast C1.
951 : /** The contrast is Strehl-normalized here.
952 : *
953 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
954 : *
955 : * \tparam imageT is an Eigen-like image type.
956 : */
957 : template <typename imageT>
958 : void C1Map( imageT &map, ///< [in] the map image to be filled in with contrast
959 : bool normStrehl );
960 :
961 : /// Calculate the residual variance due to measurement and time delay errors in phase/OPD at a spatial frequency.
962 : /** Used to calculate contrast \ref C2().
963 : *
964 : * \returns variance at (m,n).
965 : */
966 : realT C2var( realT m, ///< [in] is the spatial frequency index in u
967 : realT n ///< [in] is the spatial frequency index in v
968 : );
969 :
970 : /// Calculate the contrast due to measurement and time delay errors in phase/OPD at a spatial frequency.
971 : /** Contrast C2 is just the total variance due to time delay and measurement errors,
972 : * divided by the Strehl ratio. See Guyon (2005) \cite guyon_2005, and the updated
973 : * derivation in Males \& Guyon (2018) \cite males_guyon_2018.
974 : *
975 : * \returns C2.
976 : */
977 : realT C2( realT m, ///< [in] is the spatial frequency index in u
978 : realT n, ///< [in] is the spatial frequency index in v
979 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
980 : );
981 :
982 : /// Calculate a 2D map of contrast \ref C2().
983 : /** The contrast is Strehl-normalized here.
984 : *
985 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
986 : *
987 : * \tparam imageT is an Eigen-like image type.
988 : */
989 : template <typename imageT>
990 : void C2Map( imageT &map, ///< [in] the map image to be filled in with contrast
991 : bool normStrehl );
992 :
993 : /// Calculate the residual variance due to measurement and time delay errors in amplitude at a spatial frequency.
994 : /** Used to calculate contrast \ref C3().
995 : *
996 : * \returns variance at (m,n).
997 : */
998 : realT C3var( realT m, ///< [in] is the spatial frequency index in u
999 : realT n ///< [in] is the spatial frequency index in v
1000 : );
1001 :
1002 : /// Calculate the contrast due to measurement and time delay errors in amplitude at a spatial frequency.
1003 : /** Contrast C3 is just the total variance due to time delay and measurement errors,
1004 : * divided by the Strehl ratio. See Guyon (2005) \cite guyon_2005, and the updated
1005 : * derivation in Males \& Guyon (2018) \cite males_guyon_2018.
1006 : *
1007 : * \returns C3.
1008 : */
1009 : realT C3( realT m, ///< [in] is the spatial frequency index in u
1010 : realT n, ///< [in] is the spatial frequency index in v
1011 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
1012 : );
1013 :
1014 : /// Calculate a 2D map of contrast C3.
1015 : /** The contrast is Strehl-normalized here.
1016 : *
1017 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1018 : *
1019 : * \tparam imageT is an Eigen-like image type.
1020 : */
1021 : template <typename imageT>
1022 : void C3Map( imageT &map, ///< [in] the map image to be filled in with contrast
1023 : bool normStrehl );
1024 :
1025 : /// Calculate the residual variance due to scintilation-OPD chromaticity.
1026 : /** Used to calculate contrast \ref C4().
1027 : *
1028 : * \returns variance at (m,n).
1029 : */
1030 : realT C4var( realT m, ///< [in] is the spatial frequency index in u
1031 : realT n ///< [in] is the spatial frequency index in v
1032 : );
1033 :
1034 : /// Calculate the contrast due to scintilation-OPD chromaticity.
1035 : /** Contrast C4 is due to the chromaticity of scintillation, causing the ODP measurement to be slightly incorrect at
1036 : * the science wavelength. See Guyon (2005) \cite guyon_2005, and the updated derivation in Males \& Guyon (2018)
1037 : * \cite males_guyon_2018.
1038 : *
1039 : * \returns C4.
1040 : */
1041 : realT C4( realT m, ///< [in] is the spatial frequency index in u
1042 : realT n, ///< [in] is the spatial frequency index in v
1043 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
1044 : );
1045 :
1046 : /// Calculate a 2D map of contrast C4
1047 : /** The contrast is Strehl-normalized here.
1048 : *
1049 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1050 : *
1051 : * \tparam imageT is an Eigen-like image type.
1052 : */
1053 : template <typename imageT>
1054 : void C4Map( imageT &map, ///< [in] the map image to be filled in with contrast
1055 : bool normStrehl );
1056 :
1057 : /// Calculate the residual variance due to to scintilation-amplitude chromaticity.
1058 : /** Used to calculate contrast \ref C5().
1059 : *
1060 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1061 : *
1062 : * \returns variance at (m,n).
1063 : */
1064 : realT C5var( realT m, ///< [in] is the spatial frequency index in u
1065 : realT n ///< [in] is the spatial frequency index in v
1066 : );
1067 :
1068 : /// Calculate the contrast due to scintilation-amplitude chromaticity.
1069 : /** Contrast C5 is due to the chromaticity of scintillation, causing the amplitude measurement to be slightly
1070 : * incorrect at the science wavelength. See Guyon (2005) \cite guyon_2005, and the updated derivation in Males \&
1071 : * Guyon (2018) \cite males_guyon_2018.
1072 : *
1073 : * \returns C4.
1074 : */
1075 : realT C5( realT m, ///< [in] is the spatial frequency index in u
1076 : realT n, ///< [in] is the spatial frequency index in v
1077 : bool normStrehl = true ///< [in] flag controls whether the contrast is normalized by Strehl ratio.
1078 : );
1079 :
1080 : /// Calculate a 2D map of contrast C5
1081 : /** The contrast is Strehl-normalized here.
1082 : *
1083 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1084 : *
1085 : * \tparam imageT is an Eigen-like image type.
1086 : */
1087 : template <typename imageT>
1088 : void C5Map( imageT &map, ///< [in] the map image to be filled in with contrast
1089 : bool normStrehl );
1090 :
1091 : /// Calculate the residual variance due to chromaticity of the index of refraction of air.
1092 : /** Used to calculate contrast \ref C6().
1093 : *
1094 : * \returns variance at (m,n).
1095 : */
1096 : realT C6var( realT m, ///< [in] is the spatial frequency index in u
1097 : realT n ///< [in] is the spatial frequency index in v
1098 : );
1099 :
1100 : /// Calculate the contrast due to chromaticity of the index of refraction of air.
1101 : /** Contrast C6 is due to the index of refraction of air being wavelength dependent, causing the ODP measurement to
1102 : * be slightly incorrect at the science wavelength. See Guyon (2005) \cite guyon_2005, and the updated derivation
1103 : * in Males \& Guyon (2018) \cite males_guyon_2018.
1104 : *
1105 : * \returns C6.
1106 : */
1107 : realT C6( realT m, ///< [in] is the spatial frequency index in u
1108 : realT n, ///< [in] is the spatial frequency index in v
1109 : bool normStrehl = true ///< flag controls whether the contrast is normalized by Strehl ratio.
1110 : );
1111 :
1112 : /// Calculate a 2D map of contrast C6
1113 : /** The contrast is Strehl-normalized here.
1114 : *
1115 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1116 : *
1117 : * \tparam imageT is an Eigen-like image type.
1118 : */
1119 : template <typename imageT>
1120 : void C6Map( imageT &map, ///< [in] the map image to be filled in with contrast
1121 : bool normStrehl );
1122 :
1123 : /// Calculate the residual variance due to dispersive anisoplanatism.
1124 : /** Used to calculate contrast \ref C7().
1125 : *
1126 : * \returns variance at (m,n).
1127 : */
1128 : realT C7var( realT m, ///< [in] is the spatial frequency index in u
1129 : realT n ///< [in] is the spatial frequency index in v
1130 : );
1131 :
1132 : /// Calculate the contrast due to dispersive anisoplanatism.
1133 : /** Contrast C7 is due to atmospheric dispersion causing light at different wavelengths to take different paths
1134 : * through atmospheric turbulence. See Fitzgerald (2017, in prep).
1135 : *
1136 : * \returns C7.
1137 : */
1138 : realT C7( realT m, ///< [in] is the spatial frequency index in u
1139 : realT n, ///< [in] is the spatial frequency index in v
1140 : bool normStrehl = true ///< flag controls whether the contrast is normalized by Strehl ratio.
1141 : );
1142 :
1143 : /// Calculate a 2D map of contrast C7
1144 : /** The contrast is Strehl-normalized here.
1145 : *
1146 : * \note this is the raw PSD, it must be convolved with the PSF for the true contrast.
1147 : *
1148 : * \tparam imageT is an Eigen-like image type.
1149 : */
1150 : template <typename imageT>
1151 : void C7Map( imageT &map, ///< [in] the map image to be filled in with contrast
1152 : bool normStrehl );
1153 :
1154 : ///@}
1155 :
1156 : /// Output current parameters to a stream
1157 : /** Outputs a formatted list of all current parameters.
1158 : *
1159 : */
1160 : iosT &dumpAOSystem( iosT &ios /**< [in] a std::ostream-like stream. */ );
1161 :
1162 : /// Setup the configurator to configure this class.
1163 : void setupConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
1164 :
1165 : /// Load the configuration of this class from a configurator.
1166 : /**
1167 : * \returns `error_t::noerror` when the configured atmosphere is valid, or its typed validation error.
1168 : */
1169 : error_t loadConfig( app::appConfigurator &config /**< [in] the app::configurator object*/ );
1170 : };
1171 :
1172 : template <typename realT, class inputSpectT, typename iosT>
1173 528 : aoSystem<realT, inputSpectT, iosT>::aoSystem()
1174 : {
1175 48 : wfsBeta<wfs<realT, iosT>>( nullptr ); // allocate a default ideal WFS.
1176 48 : }
1177 :
1178 : template <typename realT, class inputSpectT, typename iosT>
1179 48 : aoSystem<realT, inputSpectT, iosT>::~aoSystem()
1180 : {
1181 48 : if( m_wfsBeta && m_ownWfsBeta )
1182 : {
1183 48 : delete m_wfsBeta;
1184 : }
1185 48 : }
1186 :
1187 : template <typename realT, class inputSpectT, typename iosT>
1188 0 : void aoSystem<realT, inputSpectT, iosT>::loadGuyon2005()
1189 : {
1190 0 : atm.loadGuyon2005();
1191 :
1192 0 : F0( 1.75e9 * 0.25 * math::pi<realT>() * 64. * 0.18 ); // Converting to photons/sec
1193 0 : lam_wfs( 0.55e-6 );
1194 0 : lam_sci( 1.6e-6 );
1195 0 : D( 8. );
1196 0 : starMag( 5 );
1197 :
1198 : // The rest of Guyon 05 is very ideal
1199 0 : npix_wfs( std::vector<realT>( { 12868 } ) );
1200 0 : ron_wfs( std::vector<realT>( { 0.0 } ) );
1201 0 : Fbg( std::vector<realT>( { 0.0 } ) );
1202 :
1203 0 : d_min( 8.0 / 1e3 ); // Allow super fine sampling
1204 0 : minTauWFS( (realT)( 1. / 1e9 ) ); // Just be super fast.
1205 0 : tauWFS( 1. / 1e9 ); // Just be super fast
1206 :
1207 0 : m_specsChanged = true;
1208 0 : m_dminChanged = true;
1209 0 : }
1210 :
1211 : template <typename realT, class inputSpectT, typename iosT>
1212 0 : void aoSystem<realT, inputSpectT, iosT>::loadMagAOX()
1213 : {
1214 0 : atm.loadLCO();
1215 0 : F0( 4.2e10 );
1216 0 : lam_wfs( 0.791e-6 );
1217 0 : lam_sci( 0.656e-6 );
1218 :
1219 0 : npix_wfs( std::vector<realT>( { 9024 } ) );
1220 0 : ron_wfs( std::vector<realT>( { 0.57 } ) );
1221 0 : Fbg( std::vector<realT>( { 0.22 } ) );
1222 :
1223 0 : d_min( std::vector<realT>( { 6.5 / 48.0 } ) );
1224 0 : minTauWFS( std::vector<realT>( { 1. / 3622. } ) );
1225 0 : tauWFS( 1. / 3622. );
1226 :
1227 0 : D( 6.5 );
1228 :
1229 0 : m_specsChanged = true;
1230 0 : m_dminChanged = true;
1231 0 : }
1232 :
1233 : template <typename realT, class inputSpectT, typename iosT>
1234 0 : void aoSystem<realT, inputSpectT, iosT>::loadGMagAOX()
1235 : {
1236 :
1237 0 : loadMagAOX();
1238 :
1239 0 : F0( 7.6e10 * 368.0 / ( 0.25 * math::pi<realT>() * 6.5 * 6.5 * ( 1. - 0.29 * 0.29 ) ) ); // Scaled up.
1240 :
1241 0 : D( 25.4 );
1242 :
1243 0 : m_specsChanged = true;
1244 0 : m_dminChanged = true;
1245 0 : }
1246 :
1247 : template <typename realT, class inputSpectT, typename iosT>
1248 39 : void aoSystem<realT, inputSpectT, iosT>::D( realT nD )
1249 : {
1250 39 : m_D = nD;
1251 39 : psd.D( m_D );
1252 :
1253 39 : m_specsChanged = true;
1254 39 : m_dminChanged = true;
1255 39 : }
1256 :
1257 : template <typename realT, class inputSpectT, typename iosT>
1258 137485 : realT aoSystem<realT, inputSpectT, iosT>::D()
1259 : {
1260 137485 : return m_D;
1261 : }
1262 :
1263 : template <typename realT, class inputSpectT, typename iosT>
1264 1 : void aoSystem<realT, inputSpectT, iosT>::d_min( const std::vector<realT> &nd )
1265 : {
1266 1 : m_d_min.assign( nd.begin(), nd.end() );
1267 1 : m_specsChanged = true;
1268 1 : m_dminChanged = true;
1269 1 : }
1270 :
1271 : template <typename realT, class inputSpectT, typename iosT>
1272 0 : void aoSystem<realT, inputSpectT, iosT>::d_min( int idx, realT nd )
1273 : {
1274 0 : if( m_d_min.size() < idx + 1 )
1275 : {
1276 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_d_min" ) );
1277 : }
1278 :
1279 0 : m_d_min[idx] = nd;
1280 :
1281 0 : m_specsChanged = true;
1282 0 : m_dminChanged = true;
1283 0 : }
1284 :
1285 : template <typename realT, class inputSpectT, typename iosT>
1286 3 : realT aoSystem<realT, inputSpectT, iosT>::d_min( size_t idx )
1287 : {
1288 3 : if( m_d_min.size() < idx + 1 )
1289 : {
1290 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_d_min" ) );
1291 : }
1292 :
1293 3 : return m_d_min[idx];
1294 : }
1295 :
1296 : template <typename realT, class inputSpectT, typename iosT>
1297 10 : std::vector<realT> aoSystem<realT, inputSpectT, iosT>::d_min()
1298 : {
1299 10 : return m_d_min;
1300 : }
1301 :
1302 : template <typename realT, class inputSpectT, typename iosT>
1303 1 : void aoSystem<realT, inputSpectT, iosT>::optd( bool od )
1304 : {
1305 1 : m_optd = od;
1306 1 : m_specsChanged = true;
1307 1 : m_dminChanged = true;
1308 1 : }
1309 :
1310 : template <typename realT, class inputSpectT, typename iosT>
1311 12 : bool aoSystem<realT, inputSpectT, iosT>::optd()
1312 : {
1313 12 : return m_optd;
1314 : }
1315 :
1316 : template <typename realT, class inputSpectT, typename iosT>
1317 1 : void aoSystem<realT, inputSpectT, iosT>::optd_delta( realT odd )
1318 : {
1319 1 : m_optd_delta = odd;
1320 1 : m_specsChanged = true;
1321 1 : m_dminChanged = true;
1322 1 : }
1323 :
1324 : template <typename realT, class inputSpectT, typename iosT>
1325 13 : realT aoSystem<realT, inputSpectT, iosT>::optd_delta()
1326 : {
1327 13 : return m_optd_delta;
1328 : }
1329 :
1330 : template <typename realT, class inputSpectT, typename iosT>
1331 : template <typename wfsT>
1332 : void aoSystem<realT, inputSpectT, iosT>::wfsBeta( const wfsT &w )
1333 : {
1334 : if( m_wfsBeta && m_ownWfsBeta )
1335 : {
1336 : delete m_wfsBeta;
1337 : m_wfsBeta = nullptr;
1338 : }
1339 :
1340 : m_wfsBeta = (wfs<realT, iosT> *)&w;
1341 : m_ownWfsBeta = false;
1342 : }
1343 :
1344 : template <typename realT, class inputSpectT, typename iosT>
1345 : template <typename wfsT>
1346 52 : void aoSystem<realT, inputSpectT, iosT>::wfsBeta( const wfsT *w )
1347 : {
1348 52 : if( m_wfsBeta && m_ownWfsBeta )
1349 : {
1350 4 : delete m_wfsBeta;
1351 4 : m_wfsBeta = nullptr;
1352 : }
1353 :
1354 52 : if( w )
1355 : {
1356 0 : m_wfsBeta = (wfs<realT, iosT> *)w;
1357 0 : m_ownWfsBeta = false;
1358 : }
1359 : else
1360 : {
1361 52 : m_wfsBeta = new wfsT;
1362 52 : m_ownWfsBeta = true;
1363 : }
1364 52 : }
1365 :
1366 : template <typename realT, class inputSpectT, typename iosT>
1367 5 : wfs<realT, iosT> *aoSystem<realT, inputSpectT, iosT>::wfsBeta()
1368 : {
1369 5 : return m_wfsBeta;
1370 : }
1371 :
1372 : template <typename realT, class inputSpectT, typename iosT>
1373 0 : realT aoSystem<realT, inputSpectT, iosT>::beta_p( realT m, realT n )
1374 : {
1375 0 : if( m_wfsBeta == 0 )
1376 : {
1377 0 : throw( mx::exception<verboseT>( error_t::paramnotset, "The WFS is not assigned." ) );
1378 : }
1379 :
1380 0 : return m_wfsBeta->beta_p( m, n, m_D, d_opt(), atm.r_0( m_lam_sci ) );
1381 : }
1382 :
1383 : template <typename realT, class inputSpectT, typename iosT>
1384 0 : realT aoSystem<realT, inputSpectT, iosT>::beta_r( realT m, realT n )
1385 : {
1386 0 : if( m_wfsBeta == 0 )
1387 : {
1388 0 : throw( mx::exception<verboseT>( error_t::paramnotset, "The WFS is not assigned." ) );
1389 : }
1390 0 : return m_wfsBeta->beta_r( m, n, m_D, d_opt(), atm.r_0( m_lam_sci ) );
1391 : }
1392 :
1393 : template <typename realT, class inputSpectT, typename iosT>
1394 0 : void aoSystem<realT, inputSpectT, iosT>::opticalGain( realT og )
1395 : {
1396 0 : m_opticalGain = og;
1397 0 : m_specsChanged = true;
1398 0 : m_dminChanged = true;
1399 0 : }
1400 :
1401 : template <typename realT, class inputSpectT, typename iosT>
1402 0 : realT aoSystem<realT, inputSpectT, iosT>::opticalGain()
1403 : {
1404 0 : return m_opticalGain;
1405 : }
1406 :
1407 : template <typename realT, class inputSpectT, typename iosT>
1408 16 : void aoSystem<realT, inputSpectT, iosT>::lam_wfs( realT nlam )
1409 : {
1410 16 : m_lam_wfs = nlam;
1411 16 : m_specsChanged = true;
1412 16 : m_dminChanged = true;
1413 16 : }
1414 :
1415 : template <typename realT, class inputSpectT, typename iosT>
1416 13 : realT aoSystem<realT, inputSpectT, iosT>::lam_wfs()
1417 : {
1418 13 : return m_lam_wfs;
1419 : }
1420 :
1421 : template <typename realT, class inputSpectT, typename iosT>
1422 1 : void aoSystem<realT, inputSpectT, iosT>::npix_wfs( const std::vector<realT> &npix )
1423 : {
1424 1 : m_npix_wfs.resize( npix.size() );
1425 3 : for( size_t n = 0; n < npix.size(); ++n )
1426 : {
1427 2 : m_npix_wfs[n] = npix[n];
1428 : }
1429 :
1430 1 : m_specsChanged = true;
1431 1 : m_dminChanged = true;
1432 1 : }
1433 :
1434 : template <typename realT, class inputSpectT, typename iosT>
1435 0 : void aoSystem<realT, inputSpectT, iosT>::npix_wfs( int idx, realT npix )
1436 : {
1437 0 : if( m_npix_wfs.size() < idx + 1 )
1438 : {
1439 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_npix_wfs" ) );
1440 : }
1441 :
1442 0 : m_npix_wfs[idx] = npix;
1443 :
1444 0 : m_specsChanged = true;
1445 0 : m_dminChanged = true;
1446 0 : }
1447 :
1448 : template <typename realT, class inputSpectT, typename iosT>
1449 4 : realT aoSystem<realT, inputSpectT, iosT>::npix_wfs( size_t idx )
1450 : {
1451 4 : if( m_npix_wfs.size() < idx + 1 )
1452 : {
1453 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_npix_wfs" ) );
1454 : }
1455 :
1456 4 : return m_npix_wfs[idx];
1457 : }
1458 :
1459 : template <typename realT, class inputSpectT, typename iosT>
1460 13 : std::vector<realT> aoSystem<realT, inputSpectT, iosT>::npix_wfs()
1461 : {
1462 13 : return m_npix_wfs;
1463 : }
1464 :
1465 : template <typename realT, class inputSpectT, typename iosT>
1466 1 : void aoSystem<realT, inputSpectT, iosT>::ron_wfs( const std::vector<realT> &nron )
1467 : {
1468 1 : m_ron_wfs.resize( nron.size() );
1469 3 : for( size_t n = 0; n < nron.size(); ++n )
1470 : {
1471 2 : m_ron_wfs[n] = nron[n];
1472 : }
1473 :
1474 1 : m_specsChanged = true;
1475 1 : m_dminChanged = true;
1476 1 : }
1477 :
1478 : template <typename realT, class inputSpectT, typename iosT>
1479 0 : void aoSystem<realT, inputSpectT, iosT>::ron_wfs( int idx, realT nron )
1480 : {
1481 0 : if( m_ron_wfs.size() < idx + 1 )
1482 : {
1483 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_ron_wfs" ) );
1484 : }
1485 :
1486 0 : m_ron_wfs[idx] = nron;
1487 :
1488 0 : m_specsChanged = true;
1489 0 : m_dminChanged = true;
1490 0 : }
1491 :
1492 : template <typename realT, class inputSpectT, typename iosT>
1493 4 : realT aoSystem<realT, inputSpectT, iosT>::ron_wfs( size_t idx )
1494 : {
1495 4 : if( m_ron_wfs.size() < idx + 1 )
1496 : {
1497 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_ron_wfs" ) );
1498 : }
1499 :
1500 4 : return m_ron_wfs[idx];
1501 : }
1502 :
1503 : template <typename realT, class inputSpectT, typename iosT>
1504 13 : std::vector<realT> aoSystem<realT, inputSpectT, iosT>::ron_wfs()
1505 : {
1506 13 : return m_ron_wfs;
1507 : }
1508 :
1509 : template <typename realT, class inputSpectT, typename iosT>
1510 1 : void aoSystem<realT, inputSpectT, iosT>::Fbg( const std::vector<realT> &fbg )
1511 : {
1512 1 : m_Fbg.resize( fbg.size() );
1513 3 : for( size_t n = 0; n < fbg.size(); ++n )
1514 : {
1515 2 : m_Fbg[n] = fbg[n];
1516 : }
1517 :
1518 1 : m_specsChanged = true;
1519 1 : m_dminChanged = true;
1520 1 : }
1521 :
1522 : template <typename realT, class inputSpectT, typename iosT>
1523 0 : void aoSystem<realT, inputSpectT, iosT>::Fbg( int idx, realT fbg )
1524 : {
1525 0 : if( m_Fbg.size() < idx + 1 )
1526 : {
1527 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_Fbg" ) );
1528 : }
1529 :
1530 0 : m_Fbg[idx] = fbg;
1531 :
1532 0 : m_specsChanged = true;
1533 0 : m_dminChanged = true;
1534 0 : }
1535 :
1536 : template <typename realT, class inputSpectT, typename iosT>
1537 4 : realT aoSystem<realT, inputSpectT, iosT>::Fbg( size_t idx )
1538 : {
1539 4 : if( m_Fbg.size() < idx + 1 )
1540 : {
1541 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_Fbg" ) );
1542 : }
1543 :
1544 4 : return m_Fbg[idx];
1545 : }
1546 :
1547 : template <typename realT, class inputSpectT, typename iosT>
1548 13 : std::vector<realT> aoSystem<realT, inputSpectT, iosT>::Fbg()
1549 : {
1550 13 : return m_Fbg;
1551 : }
1552 :
1553 : template <typename realT, class inputSpectT, typename iosT>
1554 1 : void aoSystem<realT, inputSpectT, iosT>::minTauWFS( const std::vector<realT> &ntau )
1555 : {
1556 1 : m_minTauWFS.resize( ntau.size() );
1557 3 : for( size_t n = 0; n < ntau.size(); ++n )
1558 : {
1559 2 : m_minTauWFS[n] = ntau[n];
1560 : }
1561 :
1562 1 : m_specsChanged = true;
1563 1 : m_dminChanged = true;
1564 1 : }
1565 :
1566 : template <typename realT, class inputSpectT, typename iosT>
1567 0 : void aoSystem<realT, inputSpectT, iosT>::minTauWFS( size_t idx, realT ntau )
1568 : {
1569 0 : if( m_minTauWFS.size() < idx + 1 )
1570 : {
1571 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_ntau_wfs" ) );
1572 : }
1573 :
1574 0 : m_minTauWFS[idx] = ntau;
1575 :
1576 0 : m_specsChanged = true;
1577 0 : m_dminChanged = true;
1578 0 : }
1579 :
1580 : template <typename realT, class inputSpectT, typename iosT>
1581 4 : realT aoSystem<realT, inputSpectT, iosT>::minTauWFS( size_t idx )
1582 : {
1583 4 : if( m_minTauWFS.size() < idx + 1 )
1584 : {
1585 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "idx larger than m_ntau_wfs" ) );
1586 : }
1587 :
1588 4 : return m_minTauWFS[idx];
1589 : }
1590 :
1591 : template <typename realT, class inputSpectT, typename iosT>
1592 13 : std::vector<realT> aoSystem<realT, inputSpectT, iosT>::minTauWFS()
1593 : {
1594 13 : return m_minTauWFS;
1595 : }
1596 :
1597 : template <typename realT, class inputSpectT, typename iosT>
1598 1 : void aoSystem<realT, inputSpectT, iosT>::bin_npix( bool bnp )
1599 : {
1600 1 : if( bnp != m_bin_npix )
1601 : {
1602 1 : m_bin_npix = bnp;
1603 1 : m_specsChanged = true;
1604 1 : m_dminChanged = true;
1605 : }
1606 1 : }
1607 :
1608 : template <typename realT, class inputSpectT, typename iosT>
1609 12 : bool aoSystem<realT, inputSpectT, iosT>::bin_npix()
1610 : {
1611 12 : return m_bin_npix;
1612 : }
1613 :
1614 : template <typename realT, class inputSpectT, typename iosT>
1615 0 : int aoSystem<realT, inputSpectT, iosT>::bin_opt()
1616 : {
1617 0 : d_opt();
1618 0 : return m_bin_opt;
1619 : }
1620 :
1621 : template <typename realT, class inputSpectT, typename iosT>
1622 1 : void aoSystem<realT, inputSpectT, iosT>::tauWFS( realT ntau )
1623 : {
1624 1 : m_tauWFS = ntau;
1625 1 : m_specsChanged = true;
1626 1 : m_dminChanged = true;
1627 1 : }
1628 :
1629 : template <typename realT, class inputSpectT, typename iosT>
1630 13 : realT aoSystem<realT, inputSpectT, iosT>::tauWFS()
1631 : {
1632 13 : return m_tauWFS;
1633 : }
1634 :
1635 : template <typename realT, class inputSpectT, typename iosT>
1636 1 : void aoSystem<realT, inputSpectT, iosT>::deltaTau( realT ndel )
1637 : {
1638 1 : m_deltaTau = ndel;
1639 1 : m_specsChanged = true;
1640 1 : m_dminChanged = true;
1641 1 : }
1642 :
1643 : template <typename realT, class inputSpectT, typename iosT>
1644 13 : realT aoSystem<realT, inputSpectT, iosT>::deltaTau()
1645 : {
1646 13 : return m_deltaTau;
1647 : }
1648 :
1649 : template <typename realT, class inputSpectT, typename iosT>
1650 1 : void aoSystem<realT, inputSpectT, iosT>::optTau( bool ot )
1651 : {
1652 1 : m_optTau = ot;
1653 1 : m_specsChanged = true;
1654 1 : m_dminChanged = true;
1655 1 : }
1656 :
1657 : template <typename realT, class inputSpectT, typename iosT>
1658 12 : bool aoSystem<realT, inputSpectT, iosT>::optTau()
1659 : {
1660 12 : return m_optTau;
1661 : }
1662 :
1663 : template <typename realT, class inputSpectT, typename iosT>
1664 18 : void aoSystem<realT, inputSpectT, iosT>::lam_sci( realT nlam )
1665 : {
1666 18 : m_lam_sci = nlam;
1667 18 : m_specsChanged = true;
1668 18 : m_dminChanged = true;
1669 18 : }
1670 :
1671 : template <typename realT, class inputSpectT, typename iosT>
1672 20 : realT aoSystem<realT, inputSpectT, iosT>::lam_sci()
1673 : {
1674 20 : return m_lam_sci;
1675 : }
1676 :
1677 : template <typename realT, class inputSpectT, typename iosT>
1678 1 : void aoSystem<realT, inputSpectT, iosT>::zeta( realT nz )
1679 : {
1680 1 : m_zeta = nz;
1681 1 : m_secZeta = 1 / cos( m_zeta );
1682 :
1683 1 : m_specsChanged = true;
1684 1 : m_dminChanged = true;
1685 1 : }
1686 :
1687 : template <typename realT, class inputSpectT, typename iosT>
1688 13 : realT aoSystem<realT, inputSpectT, iosT>::zeta()
1689 : {
1690 13 : return m_zeta;
1691 : }
1692 :
1693 : template <typename realT, class inputSpectT, typename iosT>
1694 96452 : realT aoSystem<realT, inputSpectT, iosT>::secZeta()
1695 : {
1696 96452 : return m_secZeta;
1697 : }
1698 :
1699 : template <typename realT, class inputSpectT, typename iosT>
1700 1 : void aoSystem<realT, inputSpectT, iosT>::fit_mn_max( int mnm )
1701 : {
1702 1 : if( mnm < 0 )
1703 0 : mnm = 0;
1704 1 : m_fit_mn_max = mnm;
1705 1 : }
1706 :
1707 : template <typename realT, class inputSpectT, typename iosT>
1708 12 : int aoSystem<realT, inputSpectT, iosT>::fit_mn_max()
1709 : {
1710 12 : return m_fit_mn_max;
1711 : }
1712 :
1713 : template <typename realT, class inputSpectT, typename iosT>
1714 1 : void aoSystem<realT, inputSpectT, iosT>::circularLimit( bool cl )
1715 : {
1716 1 : m_circularLimit = cl;
1717 1 : m_specsChanged = true;
1718 1 : m_dminChanged = true;
1719 1 : }
1720 :
1721 : template <typename realT, class inputSpectT, typename iosT>
1722 12 : bool aoSystem<realT, inputSpectT, iosT>::circularLimit()
1723 : {
1724 12 : return m_circularLimit;
1725 : }
1726 :
1727 : template <typename realT, class inputSpectT, typename iosT>
1728 1 : void aoSystem<realT, inputSpectT, iosT>::spatialFilter_ku( realT kx )
1729 : {
1730 1 : m_spatialFilter_ku = fabs( kx );
1731 1 : m_specsChanged = true; // not sure if needed
1732 1 : m_dminChanged = true; // not sure if needed
1733 1 : }
1734 :
1735 : template <typename realT, class inputSpectT, typename iosT>
1736 50 : realT aoSystem<realT, inputSpectT, iosT>::spatialFilter_ku()
1737 : {
1738 50 : return m_spatialFilter_ku;
1739 : }
1740 :
1741 : template <typename realT, class inputSpectT, typename iosT>
1742 1 : void aoSystem<realT, inputSpectT, iosT>::spatialFilter_kv( realT ky )
1743 : {
1744 1 : m_spatialFilter_kv = fabs( ky );
1745 1 : m_specsChanged = true; // not sure if needed
1746 1 : m_dminChanged = true; // not sure if needed
1747 1 : }
1748 :
1749 : template <typename realT, class inputSpectT, typename iosT>
1750 50 : realT aoSystem<realT, inputSpectT, iosT>::spatialFilter_kv()
1751 : {
1752 50 : return m_spatialFilter_kv;
1753 : }
1754 :
1755 : template <typename realT, class inputSpectT, typename iosT>
1756 1 : void aoSystem<realT, inputSpectT, iosT>::ncp_wfe( realT nwfe )
1757 : {
1758 1 : m_ncp_wfe = nwfe;
1759 1 : m_specsChanged = true;
1760 1 : m_dminChanged = true;
1761 1 : }
1762 :
1763 : template <typename realT, class inputSpectT, typename iosT>
1764 12 : realT aoSystem<realT, inputSpectT, iosT>::ncp_wfe()
1765 : {
1766 12 : return m_ncp_wfe;
1767 : }
1768 :
1769 : template <typename realT, class inputSpectT, typename iosT>
1770 1 : void aoSystem<realT, inputSpectT, iosT>::ncp_alpha( realT alpha )
1771 : {
1772 1 : m_ncp_alpha = alpha;
1773 1 : m_specsChanged = true;
1774 1 : m_dminChanged = true;
1775 1 : }
1776 :
1777 : template <typename realT, class inputSpectT, typename iosT>
1778 12 : realT aoSystem<realT, inputSpectT, iosT>::ncp_alpha()
1779 : {
1780 12 : return m_ncp_alpha;
1781 : }
1782 :
1783 : template <typename realT, class inputSpectT, typename iosT>
1784 1 : void aoSystem<realT, inputSpectT, iosT>::F0( realT nF0 )
1785 : {
1786 1 : m_F0 = nF0;
1787 1 : m_specsChanged = true;
1788 1 : m_dminChanged = true;
1789 1 : }
1790 :
1791 : template <typename realT, class inputSpectT, typename iosT>
1792 13 : realT aoSystem<realT, inputSpectT, iosT>::F0()
1793 : {
1794 13 : return m_F0;
1795 : }
1796 :
1797 : template <typename realT, class inputSpectT, typename iosT>
1798 1 : void aoSystem<realT, inputSpectT, iosT>::starMag( realT nmag )
1799 : {
1800 1 : m_starMag = nmag;
1801 1 : m_specsChanged = true;
1802 1 : m_dminChanged = true;
1803 1 : }
1804 :
1805 : template <typename realT, class inputSpectT, typename iosT>
1806 13 : realT aoSystem<realT, inputSpectT, iosT>::starMag()
1807 : {
1808 13 : return m_starMag;
1809 : }
1810 :
1811 : template <typename realT, class inputSpectT, typename iosT>
1812 4 : realT aoSystem<realT, inputSpectT, iosT>::Fg( realT mag )
1813 : {
1814 4 : return m_F0 * pow( 10.0, -0.4 * mag );
1815 : }
1816 :
1817 : template <typename realT, class inputSpectT, typename iosT>
1818 2 : realT aoSystem<realT, inputSpectT, iosT>::Fg()
1819 : {
1820 2 : return Fg( m_starMag );
1821 : }
1822 :
1823 : template <typename realT, class inputSpectT, typename iosT>
1824 0 : const Eigen::Array<int, -1, -1> &aoSystem<realT, inputSpectT, iosT>::controlledModes()
1825 : {
1826 0 : if( m_specsChanged || m_dminChanged )
1827 0 : calcStrehl();
1828 :
1829 0 : return m_controlledModes;
1830 : }
1831 :
1832 : template <typename realT, class inputSpectT, typename iosT>
1833 0 : realT aoSystem<realT, inputSpectT, iosT>::signal2Noise2( realT &Nph, realT tau_wfs, realT d, int b )
1834 : {
1835 0 : Nph = Fg() * tau_wfs;
1836 :
1837 0 : double binfact = 1.0;
1838 0 : int binidx = 0;
1839 0 : if( m_bin_npix )
1840 : {
1841 0 : if( m_npix_wfs.size() == 1 ) // Only if "true binning", otherwise the WFS mode vectors handle it.
1842 : {
1843 0 : binfact = 1. / pow( (realT)b + 1, 2 );
1844 : }
1845 : else
1846 0 : binidx = b;
1847 : }
1848 :
1849 0 : return pow( Nph, 2 ) / ( m_npix_wfs[binidx] * binfact * ( m_Fbg[binidx] * tau_wfs + pow( m_ron_wfs[binidx], 2 ) ) );
1850 : }
1851 :
1852 : template <typename realT, class inputSpectT, typename iosT>
1853 0 : realT aoSystem<realT, inputSpectT, iosT>::measurementError( realT m, realT n, realT d, int b )
1854 : {
1855 0 : if( m == 0 and n == 0 )
1856 0 : return 0;
1857 :
1858 : realT tau_wfs;
1859 :
1860 0 : if( m_optTau )
1861 0 : tau_wfs = optimumTauWFS( m, n, d, b );
1862 : else
1863 0 : tau_wfs = m_tauWFS;
1864 :
1865 0 : if( m_wfsBeta == 0 )
1866 : {
1867 0 : throw( mx::exception<verboseT>( error_t::paramnotset, "The WFS is not assigned." ) );
1868 : }
1869 :
1870 0 : realT beta_p = m_wfsBeta->beta_p( m, n, m_D, d, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
1871 0 : realT beta_r = m_wfsBeta->beta_r( m, n, m_D, d, atm.r_0( m_lam_wfs ) ); // sqrt(m_opticalGain);
1872 0 : realT Nph = 0;
1873 0 : realT snr2 = signal2Noise2( Nph, tau_wfs, d, b );
1874 :
1875 0 : return ( pow( beta_r, 2 ) / snr2 + pow( beta_p, 2 ) / Nph ) * pow( m_lam_wfs / m_lam_sci, 2 );
1876 : }
1877 :
1878 : template <typename realT, class inputSpectT, typename iosT>
1879 0 : realT aoSystem<realT, inputSpectT, iosT>::measurementErrorTotal()
1880 : {
1881 0 : if( m_specsChanged || m_dminChanged )
1882 0 : calcStrehl();
1883 0 : return m_wfeMeasurement;
1884 : }
1885 :
1886 : template <typename realT, class inputSpectT, typename iosT>
1887 0 : realT aoSystem<realT, inputSpectT, iosT>::timeDelayError( realT m, realT n, realT d, int b )
1888 : {
1889 0 : if( m == 0 and n == 0 )
1890 0 : return 0;
1891 :
1892 0 : realT k = sqrt( m * m + n * n ) / m_D;
1893 :
1894 : realT tau_wfs;
1895 :
1896 0 : if( m_optTau )
1897 0 : tau_wfs = optimumTauWFS( m, n, d, b );
1898 : else
1899 0 : tau_wfs = m_tauWFS;
1900 :
1901 0 : realT tau = tau_wfs + m_deltaTau;
1902 :
1903 : ///\todo handle multiple layers
1904 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
1905 0 : sqrt( atm.X( k, m_lam_sci, m_secZeta ) ) * pow( math::two_pi<realT>() * atm.v_wind() * k, 2 ) *
1906 0 : pow( tau, 2 );
1907 : }
1908 :
1909 : template <typename realT, class inputSpectT, typename iosT>
1910 0 : realT aoSystem<realT, inputSpectT, iosT>::timeDelayErrorTotal()
1911 : {
1912 0 : if( m_specsChanged || m_dminChanged )
1913 0 : calcStrehl();
1914 0 : return m_wfeTimeDelay;
1915 : }
1916 :
1917 : template <typename realT, class inputSpectT, typename iosT>
1918 0 : realT aoSystem<realT, inputSpectT, iosT>::fittingError( realT m, realT n )
1919 : {
1920 0 : realT k = sqrt( m * m + n * n ) / m_D;
1921 :
1922 : ///\todo handle multiple layers
1923 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 );
1924 : }
1925 :
1926 : template <typename realT, class inputSpectT, typename iosT>
1927 0 : realT aoSystem<realT, inputSpectT, iosT>::fittingErrorTotal()
1928 : {
1929 0 : if( m_specsChanged || m_dminChanged )
1930 0 : calcStrehl();
1931 0 : return m_wfeFitting;
1932 : }
1933 :
1934 : template <typename realT, class inputSpectT, typename iosT>
1935 0 : realT aoSystem<realT, inputSpectT, iosT>::chromScintOPDError( int m, int n )
1936 : {
1937 0 : return C4var( m, n );
1938 : }
1939 :
1940 : template <typename realT, class inputSpectT, typename iosT>
1941 0 : realT aoSystem<realT, inputSpectT, iosT>::chromScintOPDErrorTotal()
1942 : {
1943 0 : if( m_specsChanged || m_dminChanged )
1944 0 : calcStrehl();
1945 0 : return m_wfeChromScintOPD;
1946 : }
1947 :
1948 : template <typename realT, class inputSpectT, typename iosT>
1949 0 : realT aoSystem<realT, inputSpectT, iosT>::chromScintAmpError( int m, int n )
1950 : {
1951 0 : return C5var( m, n );
1952 :
1953 : #if 0
1954 : int mn_max = floor(0.5*m_D/d_opt());
1955 :
1956 : realT sum = 0;
1957 :
1958 : for(int m = -mn_max; m <= mn_max; ++m)
1959 : {
1960 : for(int n = -mn_max; n <= mn_max; ++n)
1961 : {
1962 : if(n == 0 && m == 0) continue;
1963 :
1964 : sum += C5var(m,n);
1965 : }
1966 : }
1967 :
1968 : return sum;
1969 : #endif
1970 : }
1971 :
1972 : template <typename realT, class inputSpectT, typename iosT>
1973 0 : realT aoSystem<realT, inputSpectT, iosT>::chromIndexError( int m, int n )
1974 : {
1975 0 : return C6var( m, n );
1976 : }
1977 :
1978 : template <typename realT, class inputSpectT, typename iosT>
1979 0 : realT aoSystem<realT, inputSpectT, iosT>::chromIndexErrorTotal()
1980 : {
1981 0 : if( m_specsChanged || m_dminChanged )
1982 0 : calcStrehl();
1983 0 : return m_wfeChromIndex;
1984 : }
1985 :
1986 : template <typename realT, class inputSpectT, typename iosT>
1987 0 : realT aoSystem<realT, inputSpectT, iosT>::dispAnisoOPDError( int m, int n )
1988 : {
1989 0 : return C7var( m, n );
1990 : }
1991 :
1992 : template <typename realT, class inputSpectT, typename iosT>
1993 0 : realT aoSystem<realT, inputSpectT, iosT>::dispAnisoOPDErrorTotal()
1994 : {
1995 0 : if( m_specsChanged || m_dminChanged )
1996 0 : calcStrehl();
1997 0 : return m_wfeAnisoOPD;
1998 : }
1999 :
2000 : template <typename realT, class inputSpectT, typename iosT>
2001 0 : realT aoSystem<realT, inputSpectT, iosT>::dispAnisoAmpError()
2002 : {
2003 0 : return 0;
2004 :
2005 : #if 0
2006 : int mn_max = floor(0.5*m_D/d_opt());
2007 :
2008 : realT sum = 0;
2009 :
2010 : for(int m = -mn_max; m <= mn_max; ++m)
2011 : {
2012 : for(int n = -mn_max; n <= mn_max; ++n)
2013 : {
2014 : if(n == 0 && m == 0) continue;
2015 :
2016 : if( m_circularLimit )
2017 : {
2018 : if( m*m + n*n > mn_max*mn_max ) continue;
2019 : }
2020 :
2021 : sum += C8var(m,n);
2022 : }
2023 : }
2024 :
2025 : return sum;
2026 : #endif
2027 : }
2028 :
2029 : template <typename realT, class inputSpectT, typename iosT>
2030 0 : realT aoSystem<realT, inputSpectT, iosT>::optimumTauWFS(
2031 : realT m,
2032 : realT n,
2033 : realT dact, // here called dact due to parameter collision with root-finding
2034 : int bbin )
2035 : {
2036 0 : if( m_D == 0 )
2037 : {
2038 0 : internal::mxlib_error_report( error_t::paramnotset, "Diameter (D) not set." );
2039 0 : return -1;
2040 : }
2041 :
2042 0 : if( m_F0 == 0 )
2043 : {
2044 0 : internal::mxlib_error_report( error_t::paramnotset, "0-mag photon flux (F0) not set." );
2045 0 : return -1;
2046 : }
2047 :
2048 0 : double binfact = 1.0;
2049 0 : int binidx = 0; // index into WFS configurations
2050 0 : if( m_bin_npix )
2051 : {
2052 0 : if( m_npix_wfs.size() > 1 )
2053 : {
2054 0 : binidx = bbin;
2055 : }
2056 : else
2057 : {
2058 0 : binfact = 1.0 / pow( (realT)( bbin + 1 ), 2 );
2059 : }
2060 : }
2061 :
2062 0 : realT k = sqrt( m * m + n * n ) / m_D;
2063 :
2064 0 : realT F = Fg();
2065 :
2066 0 : if( m_wfsBeta == 0 )
2067 : {
2068 0 : throw( mx::exception<verboseT>( error_t::paramnotset, "The WFS is not assigned." ) );
2069 : }
2070 :
2071 0 : realT beta_p = m_wfsBeta->beta_p( m, n, m_D, dact, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
2072 0 : realT beta_r = m_wfsBeta->beta_r( m, n, m_D, dact, atm.r_0( m_lam_wfs ) ) / sqrt( m_opticalGain );
2073 :
2074 : // Set up for root finding:
2075 : realT a, b, c, d, e;
2076 :
2077 : ///\todo handle multiple layers
2078 0 : realT Atmp = 2 * pow( atm.lam_0(), 2 ) * psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) *
2079 0 : ( atm.X( k, m_lam_wfs, m_secZeta ) ) * pow( math::two_pi<realT>() * atm.v_wind() * k, 2 );
2080 :
2081 0 : a = Atmp;
2082 0 : b = Atmp * m_deltaTau;
2083 0 : c = 0;
2084 0 : d = -1 * ( pow( m_lam_wfs, 2 ) / F ) *
2085 0 : ( ( m_npix_wfs[binidx] * binfact * m_Fbg[binidx] / F ) * pow( beta_r, 2 ) + pow( beta_p, 2 ) );
2086 0 : e = -2 * pow( m_lam_wfs, 2 ) * ( m_npix_wfs[binidx] * binfact ) * pow( m_ron_wfs[binidx], 2 ) * pow( beta_r, 2 ) /
2087 0 : pow( F, 2 );
2088 :
2089 0 : std::vector<std::complex<realT>> x;
2090 :
2091 : // Get the roots
2092 0 : math::quarticRoots( x, a, b, c, d, e );
2093 :
2094 : // Now pick the largest positive real root
2095 0 : realT tauopt = 0.0;
2096 :
2097 0 : for( int i = 0; i < 4; i++ )
2098 : {
2099 0 : if( real( x[i] ) > 0 && imag( x[i] ) == 0 && real( x[i] ) > tauopt )
2100 0 : tauopt = real( x[i] );
2101 : }
2102 :
2103 0 : if( tauopt < m_minTauWFS[binidx] )
2104 0 : tauopt = m_minTauWFS[binidx];
2105 :
2106 0 : return tauopt;
2107 0 : }
2108 :
2109 : template <typename realT, class inputSpectT, typename iosT>
2110 0 : realT aoSystem<realT, inputSpectT, iosT>::optimumTauWFS( realT m, realT n )
2111 : {
2112 0 : d_opt(); // also gets m_bin_opt;
2113 0 : return optimumTauWFS( m, n, m_d_opt, m_bin_opt );
2114 : }
2115 :
2116 : template <typename realT, class inputSpectT, typename iosT>
2117 1 : realT aoSystem<realT, inputSpectT, iosT>::d_opt()
2118 : {
2119 1 : if( !m_dminChanged )
2120 0 : return m_d_opt;
2121 :
2122 : // Validate the modes
2123 : ///\todo this should be in a separate valideModes function
2124 1 : if( m_npix_wfs.size() < 1 )
2125 : {
2126 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "npix_wfs must have at least one entry" ) );
2127 : }
2128 :
2129 1 : if( m_d_min.size() != m_npix_wfs.size() )
2130 : {
2131 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "d_min must be the same size as npix_wfs" ) );
2132 : }
2133 :
2134 1 : if( m_ron_wfs.size() != m_npix_wfs.size() )
2135 : {
2136 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "ron_wfs must be the same size as npix_wfs" ) );
2137 : }
2138 :
2139 1 : if( m_Fbg.size() != m_npix_wfs.size() )
2140 : {
2141 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "F_bg must be the same size as npix_wfs" ) );
2142 : }
2143 :
2144 1 : if( m_minTauWFS.size() != m_npix_wfs.size() )
2145 : {
2146 0 : throw( mx::exception<verboseT>( error_t::sizeerr, "minTauWFS must be the same size as npix_wfs" ) );
2147 : }
2148 :
2149 : ///\todo investigate why we tolerate m_d_min = 0 here. Is this just a config error?
2150 1 : if( m_d_min.size() == 1 && m_d_min[0] == 0 )
2151 : {
2152 1 : m_d_opt = 1e-50;
2153 1 : m_bin_opt = 0;
2154 1 : m_dminChanged = false;
2155 1 : return m_d_opt;
2156 : }
2157 :
2158 0 : if( !m_optd )
2159 : {
2160 0 : m_d_opt = m_d_min[0];
2161 0 : m_bin_opt = 0;
2162 0 : m_dminChanged = false;
2163 :
2164 0 : return m_d_opt;
2165 : }
2166 :
2167 0 : realT best_d = 0;
2168 :
2169 0 : if( m_bin_npix )
2170 : {
2171 0 : if( m_npix_wfs.size() > 1 ) // Optimize over the WFS modes
2172 : {
2173 0 : int best_idx = 0;
2174 0 : best_d = m_d_min[0];
2175 :
2176 0 : realT bestStrehlOverall = 0;
2177 :
2178 0 : for( int b = 0; b < m_npix_wfs.size(); ++b )
2179 : {
2180 0 : realT d = m_d_min[b];
2181 0 : realT s = strehl( d, b );
2182 0 : realT bestStrehl = s;
2183 :
2184 : // Find the actuator pitch at which AO error is less than fitting error at Nyquist
2185 0 : while( d < m_D / 2 )
2186 : {
2187 0 : d += m_d_min[b] / m_optd_delta;
2188 0 : s = strehl( d, b );
2189 :
2190 0 : if( s < bestStrehl ) // Strehl got worse. Can't be <= otherwise small m_optd_deltas will
2191 : // immediately out
2192 : {
2193 0 : d -= m_d_min[b] / m_optd_delta; // go back to previous d
2194 0 : break;
2195 : }
2196 : else
2197 0 : bestStrehl = s;
2198 : }
2199 :
2200 : // Check if this is optimum so far
2201 0 : if( bestStrehl >= bestStrehlOverall ) // >= to favor fewer actuators
2202 : {
2203 0 : bestStrehlOverall = bestStrehl;
2204 0 : best_idx = b;
2205 0 : best_d = d;
2206 : }
2207 : }
2208 :
2209 0 : m_bin_opt = best_idx;
2210 : }
2211 : else
2212 : {
2213 0 : realT d = m_d_min[0];
2214 0 : m_bin_opt = 0;
2215 :
2216 0 : realT s = strehl( d, m_bin_opt );
2217 0 : realT bestStrehl = s;
2218 :
2219 0 : while( d < m_D / 2 )
2220 : {
2221 0 : d += m_d_min[0] / m_optd_delta;
2222 0 : m_bin_opt = d / m_d_min[0] - 1;
2223 0 : if( m_bin_opt < 0 )
2224 0 : m_bin_opt = 0;
2225 :
2226 0 : s = strehl( d, m_bin_opt );
2227 :
2228 0 : if( s < bestStrehl ) // Strehl got worse. Can't be <= otherwise small m_optd_deltas will immediately
2229 : // out
2230 : {
2231 0 : d -= m_d_min[0] / m_optd_delta; // go back to previous d
2232 0 : m_bin_opt = d / m_d_min[0] - 1;
2233 0 : if( m_bin_opt < 0 )
2234 0 : m_bin_opt = 0;
2235 0 : break;
2236 : }
2237 : else
2238 0 : bestStrehl = s;
2239 : }
2240 0 : best_d = d;
2241 : }
2242 : }
2243 : else
2244 : {
2245 0 : realT d = m_d_min[0];
2246 0 : m_bin_opt = 0;
2247 :
2248 0 : realT s = strehl( d, m_bin_opt );
2249 0 : realT bestStrehl = s;
2250 :
2251 : // Find the actuator pitch which minimizes total error (AO error is less than fitting error at Nyquist)
2252 0 : while( d < m_D / 2 )
2253 : {
2254 0 : d += m_d_min[0] / m_optd_delta;
2255 :
2256 0 : s = strehl( d, m_bin_opt );
2257 0 : if( s < bestStrehl ) // Strehl got worse. Can't be <= otherwise small m_optd_deltas will immediately out
2258 : {
2259 0 : d -= m_d_min[0] / m_optd_delta; // go back to previous d
2260 0 : break;
2261 : }
2262 : else
2263 0 : bestStrehl = s;
2264 : }
2265 0 : best_d = d;
2266 : }
2267 :
2268 0 : m_d_opt = best_d;
2269 0 : m_dminChanged = false;
2270 :
2271 0 : return m_d_opt;
2272 : }
2273 :
2274 : template <typename realT, class inputSpectT, typename iosT>
2275 0 : realT aoSystem<realT, inputSpectT, iosT>::ncpError( int m, int n )
2276 : {
2277 0 : if( m == 0 and n == 0 )
2278 0 : return 0;
2279 :
2280 0 : realT k = sqrt( m * m + n * n ) / m_D;
2281 :
2282 0 : realT kmax = m_fit_mn_max / m_D;
2283 0 : realT kmin = 1. / m_D;
2284 :
2285 : realT beta;
2286 0 : if( m_ncp_alpha != 2 )
2287 : {
2288 0 : beta = ( m_ncp_alpha - 2 ) / ( math::two_pi<realT>() ) * 1. /
2289 0 : ( pow( kmin, -m_ncp_alpha + 2 ) - pow( kmax, -m_ncp_alpha + 2 ) );
2290 : }
2291 : else
2292 : {
2293 0 : beta = 1. / ( math::two_pi<realT>() * log( kmax / kmin ) );
2294 : }
2295 0 : return beta * ncpError() * pow( k, -m_ncp_alpha ) * pow( kmin, 2 );
2296 : }
2297 :
2298 : template <typename realT, class inputSpectT, typename iosT>
2299 0 : realT aoSystem<realT, inputSpectT, iosT>::ncpError()
2300 : {
2301 0 : return pow( m_ncp_wfe, 2 ) * pow( math::two_pi<realT>() / m_lam_sci, 2 );
2302 : }
2303 :
2304 : template <typename realT, class inputSpectT, typename iosT>
2305 0 : realT aoSystem<realT, inputSpectT, iosT>::strehl( realT d, int b )
2306 : {
2307 :
2308 0 : int mn_max = floor( 0.5 * m_D / d );
2309 :
2310 0 : realT wfeVar = 0;
2311 :
2312 : bool controlled;
2313 0 : for( int m = -m_fit_mn_max; m <= m_fit_mn_max; ++m )
2314 : {
2315 0 : for( int n = -m_fit_mn_max; n <= m_fit_mn_max; ++n )
2316 : {
2317 0 : controlled = false;
2318 0 : if( m_circularLimit )
2319 : {
2320 0 : if( m * m + n * n <= mn_max * mn_max )
2321 0 : controlled = true;
2322 : }
2323 : else
2324 : {
2325 0 : if( abs( m ) <= mn_max && abs( n ) <= mn_max )
2326 0 : controlled = true;
2327 : }
2328 :
2329 0 : if( controlled )
2330 : {
2331 0 : realT wfeMeasurement = measurementError( m, n, d, b );
2332 0 : realT wfeTimeDelay = timeDelayError( m, n, d, b );
2333 0 : realT wfeChromScintOPD = chromScintOPDError( m, n );
2334 0 : realT wfeChromIndex = chromIndexError( m, n );
2335 0 : realT wfeAnisoOPD = dispAnisoOPDError( m, n );
2336 :
2337 0 : realT wfeFitting = fittingError( m, n );
2338 :
2339 0 : if( wfeFitting < wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD )
2340 : {
2341 0 : wfeVar += wfeFitting;
2342 : }
2343 : else
2344 : {
2345 0 : wfeVar += wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD;
2346 : }
2347 : }
2348 : else
2349 : {
2350 0 : wfeVar += fittingError( m, n );
2351 : }
2352 : }
2353 : }
2354 :
2355 0 : wfeVar += ncpError();
2356 0 : return exp( -1 * wfeVar );
2357 : }
2358 :
2359 : template <typename realT, class inputSpectT, typename iosT>
2360 0 : void aoSystem<realT, inputSpectT, iosT>::calcStrehl()
2361 : {
2362 0 : realT d = d_opt(); // have to run this to get m_bin_opt set.
2363 0 : int b = m_bin_opt;
2364 :
2365 0 : int mn_max = floor( 0.5 * m_D / d );
2366 :
2367 0 : m_wfeMeasurement = 0;
2368 0 : m_wfeTimeDelay = 0;
2369 0 : m_wfeChromScintOPD = 0;
2370 0 : m_wfeChromIndex = 0;
2371 0 : m_wfeAnisoOPD = 0;
2372 :
2373 0 : m_wfeFitting = 0;
2374 :
2375 0 : m_wfeNCP = 0;
2376 :
2377 0 : m_controlledModes.resize( 2 * m_fit_mn_max + 1, 2 * m_fit_mn_max + 1 );
2378 0 : m_controlledModes.setZero();
2379 :
2380 : bool controlled;
2381 0 : for( int m = -m_fit_mn_max; m <= m_fit_mn_max; ++m )
2382 : {
2383 0 : for( int n = -m_fit_mn_max; n <= m_fit_mn_max; ++n )
2384 : {
2385 0 : if( m == 0 && n == 0 )
2386 : {
2387 0 : continue;
2388 : }
2389 :
2390 0 : controlled = false;
2391 0 : if( m_circularLimit )
2392 : {
2393 0 : if( m * m + n * n <= mn_max * mn_max )
2394 0 : controlled = true;
2395 : }
2396 : else
2397 : {
2398 0 : if( abs( m ) <= mn_max && abs( n ) <= mn_max )
2399 0 : controlled = true;
2400 : }
2401 :
2402 0 : if( controlled )
2403 : {
2404 0 : realT wfeMeasurement = measurementError( m, n, d, b );
2405 0 : realT wfeTimeDelay = timeDelayError( m, n, d, b );
2406 0 : realT wfeChromScintOPD = chromScintOPDError( m, n );
2407 0 : realT wfeChromIndex = chromIndexError( m, n );
2408 0 : realT wfeAnisoOPD = dispAnisoOPDError( m, n );
2409 :
2410 0 : realT wfeFitting = fittingError( m, n );
2411 :
2412 0 : if( wfeFitting < wfeMeasurement + wfeTimeDelay + wfeChromScintOPD + wfeChromIndex + wfeAnisoOPD )
2413 : {
2414 0 : m_wfeFitting += wfeFitting;
2415 : }
2416 : else // it's worth controlling this mode.
2417 : {
2418 0 : m_wfeMeasurement += wfeMeasurement;
2419 0 : m_wfeTimeDelay += wfeTimeDelay;
2420 0 : m_wfeChromScintOPD += wfeChromScintOPD;
2421 0 : m_wfeChromIndex += wfeChromIndex;
2422 0 : m_wfeAnisoOPD += wfeAnisoOPD;
2423 0 : m_controlledModes( m_fit_mn_max + m, m_fit_mn_max + n ) = 1;
2424 : }
2425 : }
2426 : else
2427 : {
2428 0 : m_wfeFitting += fittingError( m, n );
2429 : }
2430 : }
2431 : }
2432 :
2433 0 : m_wfeNCP = ncpError();
2434 :
2435 0 : m_wfeVar = m_wfeMeasurement + m_wfeTimeDelay + m_wfeFitting + m_wfeChromScintOPD + m_wfeChromIndex + m_wfeAnisoOPD +
2436 0 : m_wfeNCP;
2437 :
2438 0 : m_strehl = exp( -1 * m_wfeVar );
2439 :
2440 0 : m_specsChanged = false;
2441 0 : }
2442 :
2443 : template <typename realT, class inputSpectT, typename iosT>
2444 0 : realT aoSystem<realT, inputSpectT, iosT>::wfeVar()
2445 : {
2446 0 : if( m_specsChanged || m_dminChanged )
2447 0 : calcStrehl();
2448 :
2449 0 : return m_wfeVar;
2450 : }
2451 :
2452 : template <typename realT, class inputSpectT, typename iosT>
2453 0 : realT aoSystem<realT, inputSpectT, iosT>::strehl()
2454 : {
2455 0 : if( m_specsChanged || m_dminChanged )
2456 0 : calcStrehl();
2457 :
2458 0 : return m_strehl;
2459 : }
2460 :
2461 : template <typename realT, class inputSpectT, typename iosT>
2462 : template <typename varFuncT>
2463 0 : realT aoSystem<realT, inputSpectT, iosT>::C_( int m, int n, bool normStrehl, varFuncT varFunc, int doFittingError )
2464 : {
2465 0 : if( m == 0 && n == 0 )
2466 0 : return 0;
2467 :
2468 : // we run this to optimize regardless of whether we use it
2469 : // if already optimized then this is a minor hit
2470 0 : realT S = strehl();
2471 :
2472 0 : if( !normStrehl )
2473 0 : S = 1;
2474 :
2475 0 : if( doFittingError != FITTING_ERROR_NO )
2476 : {
2477 0 : int mn_max = m_D / ( 2. * d_opt() );
2478 :
2479 0 : if( m_controlledModes( m_fit_mn_max + m, m_fit_mn_max + n ) == false )
2480 : {
2481 0 : if( doFittingError == FITTING_ERROR_ZERO )
2482 0 : return 0;
2483 :
2484 0 : realT fe = fittingError( m, n );
2485 :
2486 0 : realT k = sqrt( m * m + n * n ) / m_D;
2487 :
2488 0 : if( doFittingError == FITTING_ERROR_X )
2489 0 : fe *= ( atm.X( k, m_lam_sci, m_secZeta ) );
2490 0 : else if( doFittingError == FITTING_ERROR_Y )
2491 0 : fe *= ( atm.Y( k, m_lam_sci, m_secZeta ) );
2492 : else
2493 : {
2494 0 : std::cerr << "Unknown doFittingError\n";
2495 0 : exit( -1 );
2496 : }
2497 :
2498 0 : return fe / S;
2499 : }
2500 : }
2501 : // get if not doing fitting error or if inside control region:
2502 :
2503 0 : realT var = ( this->*varFunc )( m, n );
2504 :
2505 0 : return var / S;
2506 : }
2507 :
2508 : template <typename realT, class inputSpectT, typename iosT>
2509 : template <typename imageT, typename CfuncT>
2510 : void aoSystem<realT, inputSpectT, iosT>::C_Map( imageT &im, CfuncT Cfunc, bool normStrehl )
2511 : {
2512 : int dim1 = im.rows();
2513 : int dim2 = im.cols();
2514 :
2515 : int mc = 0.5 * ( dim1 - 1 );
2516 : int nc = 0.5 * ( dim2 - 1 );
2517 :
2518 : for( int i = 0; i < dim1; ++i )
2519 : {
2520 : int m = i - mc;
2521 :
2522 : for( int j = 0; j < dim2; ++j )
2523 : {
2524 : int n = j - nc;
2525 :
2526 : im( i, j ) = ( this->*Cfunc )( m, n, normStrehl );
2527 : }
2528 : }
2529 : }
2530 :
2531 : template <typename realT, class inputSpectT, typename iosT>
2532 0 : realT aoSystem<realT, inputSpectT, iosT>::C0var( realT m, realT n )
2533 : {
2534 0 : realT k = sqrt( m * m + n * n ) / m_D;
2535 : ///\todo handle multiple layers
2536 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
2537 0 : ( atm.X( k, m_lam_sci, m_secZeta ) );
2538 : }
2539 :
2540 : template <typename realT, class inputSpectT, typename iosT>
2541 0 : realT aoSystem<realT, inputSpectT, iosT>::C0( realT m, realT n, bool normStrehl )
2542 : {
2543 :
2544 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C0var, FITTING_ERROR_NO );
2545 : }
2546 :
2547 : template <typename realT, class inputSpectT, typename iosT>
2548 : template <typename imageT>
2549 : void aoSystem<realT, inputSpectT, iosT>::C0Map( imageT &im, bool normStrehl )
2550 : {
2551 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C0, normStrehl );
2552 : }
2553 :
2554 : template <typename realT, class inputSpectT, typename iosT>
2555 0 : realT aoSystem<realT, inputSpectT, iosT>::C1var( realT m, realT n )
2556 : {
2557 0 : realT k = sqrt( m * m + n * n ) / m_D;
2558 :
2559 : ///\todo handle multiple layers
2560 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
2561 0 : ( atm.Y( k, m_lam_sci, m_secZeta ) );
2562 : }
2563 :
2564 : template <typename realT, class inputSpectT, typename iosT>
2565 0 : realT aoSystem<realT, inputSpectT, iosT>::C1( realT m, realT n, bool normStrehl )
2566 : {
2567 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C1var, FITTING_ERROR_NO );
2568 : }
2569 :
2570 : template <typename realT, class inputSpectT, typename iosT>
2571 : template <typename imageT>
2572 : void aoSystem<realT, inputSpectT, iosT>::C1Map( imageT &im, bool normStrehl )
2573 : {
2574 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C1, normStrehl );
2575 : }
2576 :
2577 : template <typename realT, class inputSpectT, typename iosT>
2578 0 : realT aoSystem<realT, inputSpectT, iosT>::C2var( realT m, realT n )
2579 : {
2580 0 : d_opt();
2581 0 : return measurementError( m, n, d_opt(), m_bin_opt ) + timeDelayError( m, n, d_opt(), m_bin_opt );
2582 : }
2583 :
2584 : template <typename realT, class inputSpectT, typename iosT>
2585 0 : realT aoSystem<realT, inputSpectT, iosT>::C2( realT m, realT n, bool normStrehl )
2586 : {
2587 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C2var, FITTING_ERROR_X );
2588 : }
2589 :
2590 : template <typename realT, class inputSpectT, typename iosT>
2591 : template <typename imageT>
2592 : void aoSystem<realT, inputSpectT, iosT>::C2Map( imageT &im, bool normStrehl )
2593 : {
2594 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C2, normStrehl );
2595 : }
2596 :
2597 : template <typename realT, class inputSpectT, typename iosT>
2598 0 : realT aoSystem<realT, inputSpectT, iosT>::C3var( realT m, realT n )
2599 : {
2600 0 : return 0; // measurementError(m, n) + timeDelayError(m,n);
2601 : }
2602 :
2603 : template <typename realT, class inputSpectT, typename iosT>
2604 0 : realT aoSystem<realT, inputSpectT, iosT>::C3( realT m, realT n, bool normStrehl )
2605 : {
2606 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C3var, FITTING_ERROR_ZERO ); // FITTING_ERROR_Y);
2607 : }
2608 :
2609 : template <typename realT, class inputSpectT, typename iosT>
2610 : template <typename imageT>
2611 : void aoSystem<realT, inputSpectT, iosT>::C3Map( imageT &im, bool normStrehl )
2612 : {
2613 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C3, normStrehl );
2614 : }
2615 :
2616 : template <typename realT, class inputSpectT, typename iosT>
2617 0 : realT aoSystem<realT, inputSpectT, iosT>::C4var( realT m, realT n )
2618 : {
2619 0 : realT k = sqrt( m * m + n * n ) / m_D;
2620 :
2621 : ///\todo handle multiple layers
2622 : // This does not need to be divided by X b/c we haven't multiplied by it, this is is C0/X.
2623 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
2624 0 : atm.dX( k, m_lam_sci, m_lam_wfs );
2625 : }
2626 :
2627 : template <typename realT, class inputSpectT, typename iosT>
2628 0 : realT aoSystem<realT, inputSpectT, iosT>::C4( realT m, realT n, bool normStrehl )
2629 : {
2630 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C4var, FITTING_ERROR_ZERO );
2631 : }
2632 :
2633 : template <typename realT, class inputSpectT, typename iosT>
2634 : template <typename imageT>
2635 : void aoSystem<realT, inputSpectT, iosT>::C4Map( imageT &im, bool normStrehl )
2636 : {
2637 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C4, normStrehl );
2638 : }
2639 :
2640 : template <typename realT, class inputSpectT, typename iosT>
2641 0 : realT aoSystem<realT, inputSpectT, iosT>::C5var( realT m, realT n )
2642 : {
2643 0 : realT k = sqrt( m * m + n * n ) / m_D;
2644 :
2645 : ///\todo handle multiple layers
2646 : // This does not need to be divided by Y b/c we haven't multiplied by it, this is is C1/Y.
2647 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
2648 0 : atm.dY( k, m_lam_sci, m_lam_wfs );
2649 : }
2650 :
2651 : template <typename realT, class inputSpectT, typename iosT>
2652 0 : realT aoSystem<realT, inputSpectT, iosT>::C5( realT m, realT n, bool normStrehl )
2653 : {
2654 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C5var, FITTING_ERROR_ZERO );
2655 : }
2656 :
2657 : template <typename realT, class inputSpectT, typename iosT>
2658 : template <typename imageT>
2659 : void aoSystem<realT, inputSpectT, iosT>::C5Map( imageT &im, bool normStrehl )
2660 : {
2661 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C5, normStrehl );
2662 : }
2663 :
2664 : template <typename realT, class inputSpectT, typename iosT>
2665 0 : realT aoSystem<realT, inputSpectT, iosT>::C6var( realT m, realT n )
2666 : {
2667 0 : realT ni = atm.n_air( m_lam_sci );
2668 0 : realT nw = atm.n_air( m_lam_wfs );
2669 :
2670 0 : return C0var( m, n ) * pow( ( ni - nw ) / ( ni - 1 ), 2 ); // Fixes error in Eqn 29
2671 : }
2672 :
2673 : template <typename realT, class inputSpectT, typename iosT>
2674 0 : realT aoSystem<realT, inputSpectT, iosT>::C6( realT m, realT n, bool normStrehl )
2675 : {
2676 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C6var, FITTING_ERROR_ZERO );
2677 : }
2678 :
2679 : template <typename realT, class inputSpectT, typename iosT>
2680 : template <typename imageT>
2681 : void aoSystem<realT, inputSpectT, iosT>::C6Map( imageT &im, bool normStrehl )
2682 : {
2683 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C6, normStrehl );
2684 : }
2685 :
2686 : template <typename realT, class inputSpectT, typename iosT>
2687 0 : realT aoSystem<realT, inputSpectT, iosT>::C7var( realT m, realT n )
2688 : {
2689 0 : realT k = sqrt( m * m + n * n ) / m_D;
2690 :
2691 : ///\todo this needs to handle multiple layers -- do we violate an assumption of L_0 isn't constant?
2692 0 : return psd( atm, 0, k, m_secZeta ) / pow( m_D, 2 ) * pow( atm.lam_0() / m_lam_sci, 2 ) *
2693 0 : atm.X_Z( k, m_lam_wfs, m_lam_sci, m_secZeta );
2694 : }
2695 :
2696 : template <typename realT, class inputSpectT, typename iosT>
2697 0 : realT aoSystem<realT, inputSpectT, iosT>::C7( realT m, realT n, bool normStrehl )
2698 : {
2699 0 : return C_( m, n, normStrehl, &aoSystem<realT, inputSpectT, iosT>::C7var, FITTING_ERROR_ZERO );
2700 : }
2701 :
2702 : template <typename realT, class inputSpectT, typename iosT>
2703 : template <typename imageT>
2704 : void aoSystem<realT, inputSpectT, iosT>::C7Map( imageT &im, bool normStrehl )
2705 : {
2706 : C_Map( im, &aoSystem<realT, inputSpectT, iosT>::C7, normStrehl );
2707 : }
2708 :
2709 : template <typename realT, class inputSpectT, typename iosT>
2710 1 : iosT &aoSystem<realT, inputSpectT, iosT>::dumpAOSystem( iosT &ios )
2711 : {
2712 1 : ios << "# AO Params:\n";
2713 1 : ios << "# D = " << D() << '\n';
2714 :
2715 1 : if( m_d_min.size() > 0 )
2716 : {
2717 1 : ios << "# d_min = " << d_min( (size_t)0 );
2718 1 : for( size_t n = 1; n < m_d_min.size(); ++n )
2719 0 : ios << ',' << d_min( n );
2720 1 : ios << '\n';
2721 : }
2722 : else
2723 0 : ios << "# d_min = null\n";
2724 :
2725 1 : ios << "# optd = " << std::boolalpha << m_optd << '\n';
2726 1 : ios << "# d_opt_delta = " << optd_delta() << '\n';
2727 1 : ios << "# lam_sci = " << lam_sci() << '\n';
2728 1 : ios << "# F0 = " << F0() << '\n';
2729 1 : ios << "# starMag = " << starMag() << '\n';
2730 1 : ios << "# lam_sci = " << lam_sci() << '\n';
2731 1 : ios << "# zeta = " << zeta() << '\n';
2732 1 : ios << "# lam_wfs = " << lam_wfs() << '\n';
2733 :
2734 1 : if( npix_wfs().size() > 0 )
2735 : {
2736 1 : ios << "# npix_wfs = " << npix_wfs( (size_t)0 );
2737 1 : for( size_t n = 1; n < npix_wfs().size(); ++n )
2738 0 : ios << ',' << npix_wfs( n );
2739 1 : ios << '\n';
2740 : }
2741 : else
2742 0 : ios << "# npix_wfs = null\n";
2743 :
2744 1 : if( ron_wfs().size() > 0 )
2745 : {
2746 1 : ios << "# ron_wfs = " << ron_wfs( (size_t)0 );
2747 1 : for( size_t n = 1; n < ron_wfs().size(); ++n )
2748 0 : ios << ',' << ron_wfs( n );
2749 1 : ios << '\n';
2750 : }
2751 : else
2752 0 : ios << "# ron_wfs = null\n";
2753 :
2754 1 : if( Fbg().size() > 0 )
2755 : {
2756 1 : ios << "# Fbg = " << Fbg( (size_t)0 );
2757 1 : for( size_t n = 1; n < Fbg().size(); ++n )
2758 0 : ios << ',' << Fbg( n );
2759 1 : ios << '\n';
2760 : }
2761 : else
2762 0 : ios << "# Fbg = null\n";
2763 :
2764 1 : if( minTauWFS().size() > 0 )
2765 : {
2766 1 : ios << "# minTauWFS = " << minTauWFS( (size_t)0 );
2767 1 : for( size_t n = 1; n < minTauWFS().size(); ++n )
2768 0 : ios << ',' << minTauWFS( n );
2769 1 : ios << '\n';
2770 : }
2771 : else
2772 0 : ios << "# minTauWFS = null\n";
2773 :
2774 1 : ios << "# bin_npix = " << std::boolalpha << m_bin_npix << '\n';
2775 1 : ios << "# tauWFS = " << tauWFS() << '\n';
2776 1 : ios << "# optTau = " << std::boolalpha << m_optTau << '\n';
2777 1 : ios << "# deltaTau = " << deltaTau() << '\n';
2778 1 : ios << "# fit_mn_max = " << m_fit_mn_max << '\n';
2779 1 : ios << "# spatialFilter_ku = " << m_spatialFilter_ku << '\n';
2780 1 : ios << "# spatialFilter_kv = " << m_spatialFilter_kv << '\n';
2781 1 : ios << "# ncp_wfe = " << m_ncp_wfe << '\n';
2782 1 : ios << "# ncp_alpha = " << m_ncp_alpha << '\n';
2783 :
2784 1 : if( m_wfsBeta == 0 )
2785 : {
2786 0 : throw( mx::exception<verboseT>( error_t::paramnotset, "The WFS is not assigned." ) );
2787 : }
2788 :
2789 1 : m_wfsBeta->dumpWFS( ios );
2790 1 : psd.dumpPSD( ios );
2791 1 : atm.dumpAtmosphere( ios );
2792 :
2793 1 : dumpGitStatus( ios );
2794 :
2795 1 : return ios;
2796 : }
2797 :
2798 : template <typename realT, class inputSpectT, typename iosT>
2799 10 : void aoSystem<realT, inputSpectT, iosT>::setupConfig( app::appConfigurator &config )
2800 : {
2801 : using namespace mx::app;
2802 :
2803 : // AO System configuration
2804 10 : config.add( "aosys.wfs",
2805 20 : "",
2806 20 : "aosys.wfs",
2807 : argType::Required,
2808 20 : "aosys",
2809 20 : "wfs",
2810 : false,
2811 20 : "string",
2812 30 : "The WFS type: idealWFS, unmodPyWFS, asympModPyWFS, shwfs, calculatedWFS" );
2813 10 : config.add( "aosys.wfs_beta_p",
2814 20 : "",
2815 20 : "aosys.wfs_beta_p",
2816 : argType::Required,
2817 20 : "aosys",
2818 20 : "wfs_beta_p",
2819 : false,
2820 20 : "string",
2821 30 : "The beta_p file path for calcualtedWFS" );
2822 10 : config.add( "aosys.wfs_beta_r",
2823 20 : "",
2824 20 : "aosys.wfs_beta_r",
2825 : argType::Required,
2826 20 : "aosys",
2827 20 : "wfs_beta_r",
2828 : false,
2829 20 : "string",
2830 30 : "The beta_r file path for calcualtedWFS" );
2831 10 : config.add( "aosys.wfs_sensitivity",
2832 20 : "",
2833 20 : "aosys.wfs_sensitivity",
2834 : argType::Required,
2835 20 : "aosys",
2836 20 : "wfs_sensitivity",
2837 : false,
2838 20 : "bool",
2839 30 : "Flag indicating that beta_p/beta_r are sensitivities (inverse) [default false]" );
2840 : config
2841 140 : .add( "aosys.D", "", "aosys.D", argType::Required, "aosys", "D", false, "real", "The telescope diameter [m]" );
2842 10 : config.add( "aosys.d_min",
2843 20 : "",
2844 20 : "aosys.d_min",
2845 : argType::Required,
2846 20 : "aosys",
2847 20 : "d_min",
2848 : false,
2849 20 : "real",
2850 30 : "The minimum actuator spacing [m]" );
2851 10 : config.add( "aosys.optd",
2852 20 : "",
2853 20 : "aosys.optd",
2854 : argType::Optional,
2855 20 : "aosys",
2856 20 : "optd",
2857 : false,
2858 20 : "bool",
2859 30 : "Whether or not the actuator spacing is optimized" );
2860 10 : config.add( "aosys.optd_delta",
2861 20 : "",
2862 20 : "aosys.optd_delta",
2863 : argType::Required,
2864 20 : "aosys",
2865 20 : "optd_delta",
2866 : false,
2867 20 : "bool",
2868 30 : "The fractional change from d_min used in optimization. Set to 1 (default) for integer binnings, > 1 "
2869 : "for finer sampling." );
2870 10 : config.add( "aosys.F0",
2871 20 : "",
2872 20 : "aosys.F0",
2873 : argType::Required,
2874 20 : "aosys",
2875 20 : "F0",
2876 : false,
2877 20 : "real",
2878 30 : "Zero-mag photon flux, [photons/sec]" );
2879 10 : config.add( "aosys.lam_wfs",
2880 20 : "",
2881 20 : "aosys.lam_wfs",
2882 : argType::Required,
2883 20 : "aosys",
2884 20 : "lam_wfs",
2885 : false,
2886 20 : "real",
2887 30 : "WFS wavelength [m]" );
2888 10 : config.add( "aosys.npix_wfs",
2889 20 : "",
2890 20 : "aosys.npix_wfs",
2891 : argType::Required,
2892 20 : "aosys",
2893 20 : "npix_wfs",
2894 : false,
2895 20 : "vector<real>",
2896 30 : "The number of pixels in the WFS" );
2897 10 : config.add( "aosys.ron_wfs",
2898 20 : "",
2899 20 : "aosys.ron_wfs",
2900 : argType::Required,
2901 20 : "aosys",
2902 20 : "ron_wfs",
2903 : false,
2904 20 : "vector<real>",
2905 30 : "WFS readout noise [photons/read]" );
2906 10 : config.add( "aosys.bin_npix",
2907 20 : "",
2908 20 : "aosys.bin_npix",
2909 : argType::Required,
2910 20 : "aosys",
2911 20 : "bin_npix",
2912 : false,
2913 20 : "bool",
2914 30 : "Whether or not WFS pixels are re-binned along with actuator spacing optimization" );
2915 10 : config.add( "aosys.Fbg",
2916 20 : "",
2917 20 : "aosys.Fbg",
2918 : argType::Required,
2919 20 : "aosys",
2920 20 : "Fbg",
2921 : false,
2922 20 : "vector<real>",
2923 30 : "Background counts, [counts/pix/sec]" );
2924 10 : config.add( "aosys.tauWFS",
2925 20 : "",
2926 20 : "aosys.tauWFS",
2927 : argType::Required,
2928 20 : "aosys",
2929 20 : "tauWFS",
2930 : false,
2931 20 : "real",
2932 30 : "WFS integration time [s]" );
2933 10 : config.add( "aosys.minTauWFS",
2934 20 : "",
2935 20 : "aosys.minTauWFS",
2936 : argType::Required,
2937 20 : "aosys",
2938 20 : "minTauWFS",
2939 : false,
2940 20 : "vector<real>",
2941 30 : "Minimum WFS integration time [s]" );
2942 10 : config.add( "aosys.deltaTau",
2943 20 : "",
2944 20 : "aosys.deltaTau",
2945 : argType::Required,
2946 20 : "aosys",
2947 20 : "deltaTau",
2948 : false,
2949 20 : "real",
2950 30 : "Loop delay [s]" );
2951 10 : config.add( "aosys.optTau",
2952 20 : "",
2953 20 : "aosys.optTau",
2954 : argType::Optional,
2955 20 : "aosys",
2956 20 : "optTau",
2957 : false,
2958 20 : "bool",
2959 30 : "Whether or not the integration time is optimized" );
2960 10 : config.add( "aosys.lam_sci",
2961 20 : "",
2962 20 : "aosys.lam_sci",
2963 : argType::Required,
2964 20 : "aosys",
2965 20 : "lam_sci",
2966 : false,
2967 20 : "real",
2968 30 : "Science wavelength [m]" );
2969 10 : config.add( "aosys.zeta",
2970 20 : "",
2971 20 : "aosys.zeta",
2972 : argType::Required,
2973 20 : "aosys",
2974 20 : "zeta",
2975 : false,
2976 20 : "real",
2977 30 : "Zenith distance [rad]" );
2978 10 : config.add( "aosys.fit_mn_max",
2979 20 : "",
2980 20 : "aosys.fit_mn_max",
2981 : argType::Required,
2982 20 : "aosys",
2983 20 : "fit_mn_max",
2984 : false,
2985 20 : "real",
2986 30 : "Maximum spatial frequency index to use for analysis" );
2987 10 : config.add( "aosys.circularLimit",
2988 20 : "",
2989 20 : "aosys.circularLimit",
2990 : argType::Optional,
2991 20 : "aosys",
2992 20 : "circularLimit",
2993 : false,
2994 20 : "bool",
2995 30 : " Flag to indicate that the spatial frequency limit is circular, not square." );
2996 10 : config.add( "aosys.spatialFilter_ku",
2997 20 : "",
2998 20 : "aosys.spatialFilter_ku",
2999 : argType::Required,
3000 20 : "aosys",
3001 20 : "spatialFilter_ku",
3002 : false,
3003 20 : "real",
3004 30 : "Spatial filter cutoff frequency in u [m^-1]" );
3005 10 : config.add( "aosys.spatialFilter_kv",
3006 20 : "",
3007 20 : "aosys.spatialFilter_kv",
3008 : argType::Required,
3009 20 : "aosys",
3010 20 : "spatialFilter_kv",
3011 : false,
3012 20 : "real",
3013 30 : "Spatial filter cutoff frequency in v [m^-1]" );
3014 10 : config.add( "aosys.ncp_wfe",
3015 20 : "",
3016 20 : "aosys.ncp_wfe",
3017 : argType::Required,
3018 20 : "aosys",
3019 20 : "ncp_wfe",
3020 : false,
3021 20 : "real",
3022 30 : "NCP WFE between 1 lambda/D and fit_mn_max [rad^2]" );
3023 10 : config.add( "aosys.ncp_alpha",
3024 20 : "",
3025 20 : "aosys.ncp_alpha",
3026 : argType::Required,
3027 20 : "aosys",
3028 20 : "ncp_alpha",
3029 : false,
3030 20 : "real",
3031 30 : "PSD index for NCP WFE" );
3032 10 : config.add( "aosys.starMag",
3033 20 : "",
3034 20 : "aosys.starMag",
3035 : argType::Required,
3036 20 : "aosys",
3037 20 : "starMag",
3038 : false,
3039 20 : "real",
3040 30 : "Star magnitude" );
3041 10 : config.add( "aosys.starMags",
3042 20 : "",
3043 20 : "aosys.starMags",
3044 : argType::Required,
3045 20 : "aosys",
3046 20 : "starMags",
3047 : false,
3048 20 : "real vector",
3049 20 : "A vector of star magnitudes" );
3050 :
3051 10 : atm.setupConfig( config );
3052 10 : psd.setupConfig( config );
3053 10 : }
3054 :
3055 : template <typename realT, class inputSpectT, typename iosT>
3056 10 : error_t aoSystem<realT, inputSpectT, iosT>::loadConfig( app::appConfigurator &config )
3057 : {
3058 : // WFS
3059 20 : if( config.isSet( "aosys.wfs" ) )
3060 : {
3061 4 : std::string wfsStr;
3062 4 : config( wfsStr, "aosys.wfs" );
3063 :
3064 4 : if( wfsStr == "ideal" )
3065 : {
3066 1 : wfsBeta<wfs<realT, iosT>>( nullptr );
3067 : }
3068 3 : else if( wfsStr == "unmodPyWFS" )
3069 : {
3070 1 : wfsBeta<pywfsUnmod<realT, iosT>>( nullptr );
3071 : }
3072 2 : else if( wfsStr == "asympModPyWFS" )
3073 : {
3074 1 : wfsBeta<pywfsModAsymptotic<realT>>( nullptr );
3075 : }
3076 1 : else if( wfsStr == "SHWFS" )
3077 : {
3078 1 : wfsBeta<shwfs<realT, iosT>>( nullptr );
3079 : }
3080 0 : else if( wfsStr == "calculatedWFS" )
3081 : {
3082 0 : wfsBeta<calculatedWFS<realT, iosT>>( nullptr );
3083 :
3084 0 : calculatedWFS<realT, iosT> *cwfs = static_cast<calculatedWFS<realT, iosT> *>( m_wfsBeta );
3085 0 : config( cwfs->m_beta_p_file, "aosys.wfs_beta_p" );
3086 0 : config( cwfs->m_beta_r_file, "aosys.wfs_beta_r" );
3087 0 : bool sens = cwfs->m_sensitivity;
3088 0 : config( sens, "aosys.wfs_sensitivity" );
3089 0 : if( config.isSet( "aosys.wfs_sensitivity" ) )
3090 : {
3091 0 : cwfs->m_sensitivity = sens;
3092 : }
3093 : }
3094 : else
3095 : {
3096 0 : throw( mx::exception<verboseT>( error_t::invalidarg, "unknown WFS " + wfsStr + " specified" ) );
3097 : }
3098 4 : }
3099 :
3100 : // We load the default value, get the config value (which will be the default if not set), and
3101 : // then if it was set, call the setter function so any side effects are captured.
3102 :
3103 : // diameter
3104 10 : realT nD = D();
3105 20 : config( nD, "aosys.D" );
3106 20 : if( config.isSet( "aosys.D" ) )
3107 1 : D( nD );
3108 :
3109 : // d_min
3110 10 : std::vector<realT> nd_min = d_min();
3111 20 : config( nd_min, "aosys.d_min" );
3112 20 : if( config.isSet( "aosys.d_min" ) )
3113 1 : d_min( nd_min );
3114 :
3115 10 : bool noptd = optd();
3116 20 : config( noptd, "aosys.optd" );
3117 20 : if( config.isSet( "aosys.optd" ) )
3118 1 : optd( noptd );
3119 :
3120 10 : realT noptd_delta = optd_delta();
3121 20 : config( noptd_delta, "aosys.optd_delta" );
3122 20 : if( config.isSet( "aosys.optd_delta" ) )
3123 1 : optd_delta( noptd_delta );
3124 :
3125 10 : realT nlam_wfs = lam_wfs();
3126 20 : config( nlam_wfs, "aosys.lam_wfs" );
3127 20 : if( config.isSet( "aosys.lam_wfs" ) )
3128 1 : lam_wfs( nlam_wfs );
3129 :
3130 : // npix_wfs
3131 10 : std::vector<realT> nnpix_wfs = npix_wfs();
3132 20 : config( nnpix_wfs, "aosys.npix_wfs" );
3133 20 : if( config.isSet( "aosys.npix_wfs" ) )
3134 1 : npix_wfs( nnpix_wfs );
3135 :
3136 : // ron_wfs
3137 10 : std::vector<realT> nron_wfs = ron_wfs();
3138 20 : config( nron_wfs, "aosys.ron_wfs" );
3139 20 : if( config.isSet( "aosys.ron_wfs" ) )
3140 1 : ron_wfs( nron_wfs );
3141 :
3142 : // Fbg
3143 10 : std::vector<realT> nFbg = Fbg();
3144 20 : config( nFbg, "aosys.Fbg" );
3145 20 : if( config.isSet( "aosys.Fbg" ) )
3146 1 : Fbg( nFbg );
3147 :
3148 : // minTauWFS
3149 10 : std::vector<realT> nminTauWFS = minTauWFS();
3150 20 : config( nminTauWFS, "aosys.minTauWFS" );
3151 20 : if( config.isSet( "aosys.minTauWFS" ) )
3152 1 : minTauWFS( nminTauWFS );
3153 :
3154 : // bin_npix
3155 10 : bool nbin_npix = bin_npix();
3156 20 : config( nbin_npix, "aosys.bin_npix" );
3157 20 : if( config.isSet( "aosys.bin_npix" ) )
3158 1 : bin_npix( nbin_npix );
3159 :
3160 : // tauWFS
3161 10 : realT ntauWFS = tauWFS();
3162 20 : config( ntauWFS, "aosys.tauWFS" );
3163 20 : if( config.isSet( "aosys.tauWFS" ) )
3164 1 : tauWFS( ntauWFS );
3165 :
3166 : // deltaTau
3167 10 : realT ndeltaTau = deltaTau();
3168 20 : config( ndeltaTau, "aosys.deltaTau" );
3169 20 : if( config.isSet( "aosys.deltaTau" ) )
3170 1 : deltaTau( ndeltaTau );
3171 :
3172 : // optTau
3173 10 : bool noptTau = optTau();
3174 20 : config( noptTau, "aosys.optTau" );
3175 20 : if( config.isSet( "aosys.optTau" ) )
3176 1 : optTau( noptTau );
3177 :
3178 : // lam_sci
3179 10 : realT nlam_sci = lam_sci();
3180 20 : config( nlam_sci, "aosys.lam_sci" );
3181 20 : if( config.isSet( "aosys.lam_sci" ) )
3182 1 : lam_sci( nlam_sci );
3183 :
3184 : // zeta
3185 10 : realT nzeta = zeta();
3186 20 : config( nzeta, "aosys.zeta" );
3187 20 : if( config.isSet( "aosys.zeta" ) )
3188 1 : zeta( nzeta );
3189 :
3190 : // fit_mn_max
3191 10 : realT fmnm = fit_mn_max();
3192 20 : config( fmnm, "aosys.fit_mn_max" );
3193 20 : if( config.isSet( "aosys.fit_mn_max" ) )
3194 1 : fit_mn_max( fmnm );
3195 :
3196 : // circularLimit
3197 10 : bool cl = circularLimit();
3198 20 : config( cl, "aosys.circularLimit" );
3199 20 : if( config.isSet( "aosys.circularLimit" ) )
3200 1 : circularLimit( cl );
3201 :
3202 : // spatialFilter_ku
3203 10 : realT ku = spatialFilter_ku();
3204 20 : config( ku, "aosys.spatialFilter_ku" );
3205 20 : if( config.isSet( "aosys.spatialFilter_ku" ) )
3206 1 : spatialFilter_ku( ku );
3207 :
3208 : // spatialFilter_kv
3209 10 : realT kv = spatialFilter_kv();
3210 20 : config( kv, "aosys.spatialFilter_kv" );
3211 20 : if( config.isSet( "aosys.spatialFilter_kv" ) )
3212 1 : spatialFilter_kv( kv );
3213 :
3214 : // ncp_wfe
3215 10 : realT nwfe = ncp_wfe();
3216 20 : config( nwfe, "aosys.ncp_wfe" );
3217 20 : if( config.isSet( "aosys.ncp_wfe" ) )
3218 1 : ncp_wfe( nwfe );
3219 :
3220 : // ncp_alpha
3221 10 : realT na = ncp_alpha();
3222 20 : config( na, "aosys.ncp_alpha" );
3223 20 : if( config.isSet( "aosys.ncp_alpha" ) )
3224 1 : ncp_alpha( na );
3225 :
3226 : // F0
3227 10 : realT nF0 = F0();
3228 20 : config( nF0, "aosys.F0" );
3229 20 : if( config.isSet( "aosys.F0" ) )
3230 1 : F0( nF0 );
3231 :
3232 : // star_mag
3233 10 : realT smag = starMag();
3234 20 : config( smag, "aosys.starMag" );
3235 20 : if( config.isSet( "aosys.starMag" ) )
3236 1 : starMag( smag );
3237 :
3238 10 : const error_t atmosphereStatus = atm.loadConfig( config );
3239 10 : psd.loadConfig( config );
3240 :
3241 10 : return atmosphereStatus;
3242 10 : }
3243 :
3244 : extern template class aoSystem<float, vonKarmanSpectrum<float>, std::ostream>;
3245 :
3246 : extern template class aoSystem<double, vonKarmanSpectrum<double>, std::ostream>;
3247 :
3248 : extern template class aoSystem<long double, vonKarmanSpectrum<long double>, std::ostream>;
3249 :
3250 : #ifdef HASQUAD
3251 : extern template class aoSystem<__float128, vonKarmanSpectrum<__float128>, std::ostream>;
3252 : #endif
3253 :
3254 : } // namespace analysis
3255 : } // namespace AO
3256 : } // namespace mx
3257 :
3258 : #endif // aoSystem_hpp
|
