mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
fourierTemporalPSD.hpp
Go to the documentation of this file.
1/** \file fourierTemporalPSD.hpp
2 * \author Jared R. Males (jaredmales@gmail.com)
3 * \brief Calculation of the temporal PSD of Fourier modes.
4 * \ingroup mxAOm_files
5 *
6 */
7
8//***********************************************************************//
9// Copyright 2016-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 fourierTemporalPSD_hpp
28#define fourierTemporalPSD_hpp
29
30#include <iostream>
31#include <fstream>
32
33#include <sys/stat.h>
34
35#include <gsl/gsl_integration.h>
36#include <gsl/gsl_errno.h>
37
38#include <Eigen/Dense>
39
40#include "../../math/constants.hpp"
41#include "../../math/func/jinc.hpp"
42#include "../../math/func/airyPattern.hpp"
43#include "../../math/vectorUtils.hpp"
44#include "../../ioutils/fits/fitsFile.hpp"
45#include "../../sigproc/fourierModes.hpp"
46#include "../../sigproc/psdVarMean.hpp"
47#include "../../ioutils/stringUtils.hpp"
48#include "../../ioutils/readColumns.hpp"
49#include "../../ioutils/binVector.hpp"
50#include "../../ioutils/fileUtils.hpp"
51
52#include "../../ipc/ompLoopWatcher.hpp"
53#include "../../mxError.hpp"
54
55#include "aoSystem.hpp"
56#include "aoPSDs.hpp"
57#include "wfsNoisePSD.hpp"
58#include "clAOLinearPredictor.hpp"
59#include "clGainOpt.hpp"
60#include "varmapToImage.hpp"
61#include "speckleAmpPSD.hpp"
62
63#include "aoConstants.hpp"
64
65namespace mx
66{
67namespace AO
68{
69namespace analysis
70{
71
72#ifndef WSZ
73
74 /** \def WFZ
75 * \brief Size of the GSL integration workspace
76 */
77 #define WSZ 100000
78
79#endif
80
81enum basis : unsigned int
82{
83 basic, ///< The basic sine and cosine Fourier modes
84 modified ///< The modified Fourier basis from \cite males_guyon_2017
85};
86
87// Forward declaration
88template <typename realT, typename aosysT>
89realT F_basic( realT kv, void *params );
90
91// Forward declaration
92template <typename realT, typename aosysT>
93realT F_mod( realT kv, void *params );
94
95
96/// Class to manage the calculation of temporal PSDs of the Fourier modes in atmospheric turbulence.
97/** Works with both basic (sines/cosines) and modified Fourier modes.
98 *
99 * \tparam realT is a real floating point type for calculations. Currently must be double due to gsl_integration.
100 * \tparam aosysT is an AO system type, usually of type ao_system.
101 *
102 * \todo Split off the integration parameters in a separate structure.
103 * \todo once integration parameters are in a separate structure, make this a class with protected members.
104 * \todo GSL error handling
105 *
106 * \ingroup mxAOAnalytic
107 */
108template <typename _realT, typename aosysT>
109struct fourierTemporalPSD
110{
111 /// The type for arithmetic
112 typedef _realT realT;
113
114 /// The complex type for arithmetic
115 typedef std::complex<realT> complexT;
116
117 /// Pointer to an AO system structure.
118 aosysT *m_aosys{ nullptr };
119
120 realT m_f{ 0 }; ///< the current temporal frequency
121 realT m_m{ 0 }; ///< the spatial frequency m index
122 realT m_n{ 0 }; ///< the spatial frequency n index
123 realT m_cq{ 0 }; ///< The cosine of the wind direction
124 realT m_sq{ 0 }; ///< The sine of the wind direction
125 realT m_spatialFilter{ false }; ///< Flag indicating if a spatial filter is applied
126
127 bool m_strehlOG{ false };
128 bool m_uncorrectedOG{ false };
129
130 realT m_f0{ 0 }; ///< the Berdja boiling parameter
131
132 int m_p{ 1 }; ///< The parity of the mode, +/- 1. If _useBasis==MXAO_FTPSD_BASIS_BASIC then +1 indicates cosine, -1
133 ///< indicates sine.
134 int _layer_i; ///< The index of the current layer.
135
136 int _useBasis; ///< Set to MXAO_FTPSD_BASIS_BASIC/MODIFIED/PROJECTED_* to use the basic sin/cos modes, the modified
137 ///< Fourier modes, or a projection of them.
138
139 /// Workspace for the gsl integrators, allocated to WSZ if constructed as worker (with allocate == true).
140 gsl_integration_workspace *_w;
141
142 realT _absTol; ///< The absolute tolerance to use in the GSL integrator
143 realT _relTol; ///< The relative tolerance to use in the GSL integrator
144
145 int m_mode_i; ///< Projected basis mode index
146
147 Eigen::Array<realT, -1, -1> m_modeCoeffs; ///< Coeeficients of the projection onto the Fourier modes
148 realT m_minCoeffVal;
149
150 std::vector<realT> Jps;
151 std::vector<realT> Jms;
152 std::vector<int> ps;
153 std::vector<realT> ms;
154 std::vector<realT> ns;
155
156 void initProjection()
157 {
158 Jps.resize( m_modeCoeffs.cols() );
159 Jms.resize( m_modeCoeffs.cols() );
160 ps.resize( m_modeCoeffs.cols() );
161 ms.resize( m_modeCoeffs.cols() );
162 ns.resize( m_modeCoeffs.cols() );
163
164 for( int i = 0; i < m_modeCoeffs.cols(); ++i )
165 {
166 int m, n, p;
167 sigproc::fourierModeCoordinates( m, n, p, i );
168 ps[i] = p;
169 ms[i] = m;
170 ns[i] = n;
171 }
172 }
173
174 public:
175 /// Default c'tor
176 fourierTemporalPSD();
177
178 /// Constructor with workspace allocation
179 /**
180 * \param allocate if true, then the workspace for GSL integrators is allocated.
181 */
182 explicit fourierTemporalPSD( bool allocate );
183
184 /// Destructor
185 /** Frees GSL workspace if it was allocated.
186 */
187 ~fourierTemporalPSD();
188
189 protected:
190 /// Initialize parameters to default values.
191 void initialize();
192
193 public:
194 /** \name GSL Integration Tolerances
195 * For good results it seems that absolute tolerance (absTol) needs to be 1e-10. Lower tolerances cause some
196 * frequencies to drop out, etc. Relative tolerance (relTol) seems to be less sensitive, and 1e-4 works on cases
197 * tested as of 1 Jan, 2017.
198 *
199 * See the documentation for the GSL Library integrators at
200 * (https://www.gnu.org/software/gsl/manual/htmlm_node/QAGI-adaptive-integration-on-infinite-intervals.html)
201 * @{
202 */
203
204 /// Set absolute tolerance
205 /**
206 * \param at is the new absolute tolerance.
207 */
208 void absTol( realT at );
209
210 /// Get the current absolute tolerance
211 /**
212 * \returns _absTol
213 */
214 realT absTol();
215
216 /// Set relative tolerance
217 /**
218 * \param rt is the new relative tolerance.
219 */
220 void relTol( realT rt );
221
222 /// Get the current relative tolerance
223 /**
224 * \returns _relTol
225 */
226 realT relTol();
227
228 ///@}
229
230 /// Determine the frequency of the highest V-dot-k peak
231 /**
232 * \param m the spatial frequency u index
233 * \param n the spatial frequency v index
234 *
235 * \return the frequency of the fastest peak
236 */
237 realT fastestPeak( int m, int n );
238
239 /// Calculate the temporal PSD for a Fourier mode for a single layer.
240 /**
241 *
242 * \todo implement error checking.
243 * \todo need a way to track convergence failures in integral without throwing an error.
244 * \todo need better handling of averaging for the -17/3 extension.
245 *
246 */
247 int
248 singleLayerPSD( std::vector<realT> &PSD, ///< [out] the calculated PSD
249 std::vector<realT> &freq, ///< [in] the populated temporal frequency grid defining the frequencies
250 ///< at which the PSD is calculated
251 realT m, ///< [in] the first index of the spatial frequency
252 realT n, ///< [in] the second index of the spatial frequency
253 int layer_i, ///< [in] the index of the layer, for accessing the atmosphere parameters
254 int p, ///< [in] sets which mode is calculated (if basic modes, p = -1 for sine, p = +1 for cosine)
255 realT fmax = 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The PSD
256 ///< is filled in with a -17/3 power law past this frequency.
257 );
258
259 ///\cond multilayerm_parallel
260 // Conditional to exclude from Doxygen.
261
262 protected:
263 // Type to allow overloading of the multiLayerPSD workers based on whether they are parallelized or not.
264 template <bool m_parallel>
265 struct isParallel
266 {
267 };
268
269 // Parallelized version of multiLayerPSD, with OMP directives.
270 int m_multiLayerPSD( std::vector<realT> &PSD,
271 std::vector<realT> &freq,
272 realT m,
273 realT n,
274 int p,
275 realT fmax,
276 isParallel<true> parallel );
277
278 // Non-Parallelized version of multiLayerPSD, without OMP directives.
279 int m_multiLayerPSD( std::vector<realT> &PSD,
280 std::vector<realT> &freq,
281 realT m,
282 realT n,
283 int p,
284 realT fmax,
285 isParallel<false> parallel );
286
287 ///\endcond
288
289 public:
290 /// Calculate the temporal PSD for a Fourier mode in a multi-layer model.
291 /**
292 *
293 * \tparam parallel controls whether layers are calculated in parallel. Default is true. Set to false if this is
294 * called inside a parallelized loop, as in \ref makePSDGrid.
295 *
296 * \todo implement error checking
297 * \todo handle reports of convergence failures form singleLayerPSD when implemented.
298 *
299 */
300 template <bool parallel = true>
301 int multiLayerPSD(
302 std::vector<realT> &PSD, ///< [out] the calculated PSD
303 std::vector<realT>
304 &freq, ///< [in] the populated temporal frequency grid defining at which frequencies the PSD is calculated
305 realT m, ///< [in] the first index of the spatial frequency
306 realT n, ///< [in] the second index of the spatial frequency
307 int p, ///< [in] sets which mode is calculated (if basic modes, p = -1 for sine, p = +1 for cosine)
308 realT fmax =
309 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The PSD is filled in
310 /// with a -17/3 power law past this frequency. If 0, then it is taken to be 150 Hz + 2*fastestPeak(m,n).
311 );
312
313 /// Calculate PSDs over a grid of spatial frequencies.
314 /** The grid of spatial frequencies is square, set by the maximum value of m and n.
315 *
316 * The PSDs are written as mx::binVector binary files to a directory. We do not use FITS since
317 * this adds overhead and cfitisio handles parallelization poorly due to the limitation on number of created file
318 * pointers.
319 *
320 */
321 void makePSDGrid( const std::string &dir, ///< [in] the directory for output of the PSDs
322 int mnMax, ///< [in] the maximum value of m and n in the grid.
323 realT dFreq, ///< [in] the temporal frequency spacing.
324 realT maxFreq, ///< [in] the maximum temporal frequency to calculate
325 realT fmax = 0 ///< [in] [optional] set the maximum temporal frequency for the calculation. The
326 ///< PSD is filled in with a -17/3 power law past
327 /// this frequency. If 0, then it is taken to be 150 Hz + 2*fastestPeak(m,n).
328 );
329
330 /// Analyze a PSD grid under closed-loop control.
331 /** This always analyzes the simple integrator, and can also analyze the linear predictor controller. Outputs maps
332 * of optimum gains, predictor coefficients, variances, and contrasts for a range of guide star magnitudes.
333 * Optionally calculates speckle lifetimes. Optionally writes the closed-loop PSDs and transfer functions.
334 */
335 int analyzePSDGrid( const std::string &subDir, /**< [out] the sub-directory of psdDir where to write the
336 results. Is created. */
337 const std::string &psdDir, ///< [in] the directory containing the grid of PSDs.
338 int mnMax, ///< [in] the maximum value of m and n in the grid.
339 int mnCon, ///< [in] the maximum value of m and n which can be controlled.
340 realT gfixed, ///< [in] if > 0 then this fixed gain is used in the SI.
341 int lpNc, /**< [in] the number of linear predictor coefficients to analyze.
342 If 0 then LP is not analyzed.*/
343 realT lpRegPrecision, /**< [in] the initial precision for the LP regularization
344 algorithm. Normal value is 2. Higher is faster. Decrease
345 if getting stuck in local minima.*/
346 std::vector<realT> &mags, ///< [in] the guide star magnitudes to analyze for.
347 int lifetimeTrials = 0, /**< [in] [optional] number of trials used for calculating speckle
348 lifetimes. If 0,lifetimes are not calculated. */
349 bool ucLifeTs = false, /**< [in] [optional] flag controlling whether lifetimes are
350 calculated for uncontrolled modes.*/
351 bool writePSDs = false, ///< [in] [optional] flag controlling if resultant PSDs are saved
352 bool writeXfer = false /**< [in] [optional] flag controlling if resultant
353 transfer functions are saved*/
354 );
355
356 int intensityPSD( const std::string &subDir, // sub-directory of psdDir which contains the controlled system
357 // results, and where the lifetimes will be written.
358 const std::string &psdDir, // directory containing the grid of PSDs
359 const std::string &CvdPath, // path to the covariance decomposition
360 int mnMax, ///< [in] the maximum value of m and n in the grid.
361 int mnCon, ///< [in] the maximum value of m and n which can be controlled.
362 std::vector<realT> &mags, ///< [in] the guide star magnitudes
363 int lifetimeTrials, /**< [in] [optional] number of trials used for calculating
364 speckle lifetimes. If 0, lifetimes are not
365 calculated.*/
366 bool writePSDs /**< [in] [optional] flag controlling if resultant
367 PSDs are saved*/ );
368
369 /** \name Disk Storage
370 * These methods handle writing to and reading from disk. The calculated PSDs are store in the mx::BinVector binary
371 * format.
372 *
373 * A grid of PSDs is specified by its directory name. The directory contains one frequency file (freq.binv), and a
374 * set of PSD files, named according to psd_<m>_<n>_.binv.
375 *
376 *
377 * @{
378 */
379 /// Get the frequency scale for a PSD grid.
380 /**
381 */
382 int getGridFreq( std::vector<realT> &freq, ///< [out] the vector to populate with the frequency scale.
383 const std::string &dir ///< [in] specifies the directory containing the grid.
384 );
385
386 /// Get a single PSD from a PSD grid.
387 /**
388 */
389 int getGridPSD( std::vector<realT> &psd, ///< [out] the vector to populate with the PSD.
390 const std::string &dir, ///< [in] specifies the directory containing the grid.
391 int m, ///< [in] specifies the u component of spatial frequency.
392 int n ///< [in] specifies the v component of spatial frequency.
393 );
394
395 /// Get both the frequency scale and a single PSD from a PSD grid.
396 /**
397 */
398 int getGridPSD( std::vector<realT> &freq, ///< [out] the vector to populate with the frequency scale.
399 std::vector<realT> &psd, ///< [out] the vector to populate with the PSD.
400 const std::string &dir, ///< [in] specifies the directory containing the grid.
401 int m, ///< [in] specifies the u component of spatial frequency.
402 int n ///< [in] specifies the v component of spatial frequency.
403 );
404
405 ///@}
406};
407
408template <typename realT, typename aosysT>
409fourierTemporalPSD<realT, aosysT>::fourierTemporalPSD()
410{
411 m_aosys = nullptr;
412 initialize();
413}
414
415template <typename realT, typename aosysT>
416fourierTemporalPSD<realT, aosysT>::fourierTemporalPSD( bool allocate )
417{
418 m_aosys = nullptr;
419 initialize();
420
421 if( allocate )
422 {
423 _w = gsl_integration_workspace_alloc( WSZ );
424 }
425}
426
427template <typename realT, typename aosysT>
428fourierTemporalPSD<realT, aosysT>::~fourierTemporalPSD()
429{
430 if( _w )
431 {
432 gsl_integration_workspace_free( _w );
433 }
434}
435
436template <typename realT, typename aosysT>
437void fourierTemporalPSD<realT, aosysT>::initialize()
438{
439 _useBasis = basis::modified;
440 _w = 0;
441
442 _absTol = 1e-10;
443 _relTol = 1e-4;
444}
445
446template <typename realT, typename aosysT>
447void fourierTemporalPSD<realT, aosysT>::absTol( realT at )
448{
449 _absTol = at;
450}
451
452template <typename realT, typename aosysT>
453realT fourierTemporalPSD<realT, aosysT>::absTol()
454{
455 return _absTol;
456}
457
458template <typename realT, typename aosysT>
459void fourierTemporalPSD<realT, aosysT>::relTol( realT rt )
460{
461 _relTol = rt;
462}
463
464template <typename realT, typename aosysT>
465realT fourierTemporalPSD<realT, aosysT>::relTol()
466{
467 return _relTol;
468}
469
470template <typename realT, typename aosysT>
471realT fourierTemporalPSD<realT, aosysT>::fastestPeak( int m, int n )
472{
473 realT ku, kv, vu, vv;
474
475 ku = ( (realT)m / m_aosys->D() );
476 kv = ( (realT)n / m_aosys->D() );
477
478 realT f, fmax = 0;
479
480 for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
481 {
482 vu = m_aosys->atm.layer_v_wind( i ) * cos( m_aosys->atm.layer_dir( i ) );
483 vv = m_aosys->atm.layer_v_wind( i ) * sin( m_aosys->atm.layer_dir( i ) );
484
485 f = fabs( ku * vu + kv * vv );
486 if( f > fmax )
487 fmax = f;
488 }
489
490 return fmax;
491}
492
493template <typename realT, typename aosysT>
494int fourierTemporalPSD<realT, aosysT>::singleLayerPSD(
495 std::vector<realT> &PSD, std::vector<realT> &freq, realT m, realT n, int layer_i, int p, realT fmax )
496{
497 if( fmax == 0 )
498 fmax = freq[freq.size() - 1];
499
500 realT v_wind = m_aosys->atm.layer_v_wind( layer_i );
501 realT q_wind = m_aosys->atm.layer_dir( layer_i );
502
503 // Rotate the basis
504 realT cq = cos( q_wind );
505 realT sq = sin( q_wind );
506
507 realT scale = 2 * ( 1 / v_wind ); // Factor of 2 for negative frequencies.
508
509 // We'll get the occasional failure to reach tolerance error, just ignore them all for now.
510 gsl_set_error_handler_off();
511
512 // Create a local instance so that we're reentrant
513 fourierTemporalPSD<realT, aosysT> params( true );
514
515 params.m_aosys = m_aosys;
516 params._layer_i = layer_i;
517 params.m_m = m * cq + n * sq;
518 params.m_n = -m * sq + n * cq;
519 params.m_cq = cq; // for de-rotating ku and kv for spatial filtering
520 params.m_sq = sq; // for de-rotation ku and kv for spatial filtering
521 if( m_aosys->spatialFilter_ku() < std::numeric_limits<realT>::max() ||
522 m_aosys->spatialFilter_kv() < std::numeric_limits<realT>::max() )
523 params.m_spatialFilter = true;
524
525 params.m_p = p;
526 params.m_f0 = m_f0;
527
528 params.m_mode_i = m_mode_i;
529 params.m_modeCoeffs = m_modeCoeffs;
530 params.m_minCoeffVal = m_minCoeffVal;
531
532 realT result, error;
533
534 // Setup the GSL calculation
535 gsl_function func;
536 switch( _useBasis )
537 {
538 case basis::basic: // MXAO_FTPSD_BASIS_BASIC:
539 func.function = &F_basic<realT, aosysT>;
540 break;
541 case basis::modified: // MXAO_FTPSD_BASIS_MODIFIED:
542 func.function = &F_mod<realT, aosysT>;
543 break;
544 default:
545 mxError( "fourierTemporalPSD::singleLayerPSD", MXE_INVALIDARG, "value of _useBasis is not valid." );
546 return -1;
547 }
548
549 func.params = &params;
550
551 // Here we only calculate up to fmax.
552 size_t i = 0;
553 while( freq[i] <= fmax )
554 {
555 params.m_f = freq[i];
556
557 int ec = gsl_integration_qagi( &func, _absTol, _relTol, WSZ, params._w, &result, &error );
558
559 if( ec == GSL_EDIVERGE )
560 {
561 std::cerr << "GSL_EDIVERGE:" << p << " " << freq[i] << " " << v_wind << " " << m << " " << n << " " << m_m
562 << " " << m_n << "\n";
563 std::cerr << "ignoring . . .\n";
564 }
565
566 PSD[i] = scale * result;
567
568 ++i;
569 if( i >= freq.size() )
570 break;
571 }
572
573 // Now fill in from fmax to the actual max frequency with a -(alpha+2) power law.
574 size_t j = i;
575
576 if( j == freq.size() )
577 return 0;
578
579 // First average result for last 50.
580 PSD[j] =
581 PSD[i - 50] * pow( freq[i - 50] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds<realT>());
582 for( size_t k = 49; k > 0; --k )
583 {
584 PSD[j] +=
585 PSD[i - k] * pow( freq[i - k] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds<realT>());
586 }
587 PSD[j] /= 50.0;
588 ++j;
589 ++i;
590 if( j == freq.size() )
591 return 0;
592 while( j < freq.size() )
593 {
594 PSD[j] =
595 PSD[i - 1] * pow( freq[i - 1] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds<realT>());
596 ++j;
597 }
598
599 return 0;
600}
601
602template <typename realT, typename aosysT>
603int fourierTemporalPSD<realT, aosysT>::m_multiLayerPSD(
604 std::vector<realT> &PSD, std::vector<realT> &freq, realT m, realT n, int p, realT fmax, isParallel<true> parallel )
605{
606 static_cast<void>( parallel );
607
608#pragma omp parallel
609 {
610 // Records each layer PSD
611 std::vector<realT> single_PSD( freq.size() );
612
613#pragma omp for
614 for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
615 {
616 singleLayerPSD( single_PSD, freq, m, n, i, p, fmax );
617
618// Now add the single layer PSD to the overall PSD, weighted by Cn2
619#pragma omp critical
620 for( size_t j = 0; j < freq.size(); ++j )
621 {
622 PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j];
623 }
624 }
625 }
626
627 return 0;
628}
629
630template <typename realT, typename aosysT>
631int fourierTemporalPSD<realT, aosysT>::m_multiLayerPSD(
632 std::vector<realT> &PSD, std::vector<realT> &freq, realT m, realT n, int p, realT fmax, isParallel<false> parallel )
633{
634 static_cast<void>( parallel );
635
636 // Records each layer PSD
637 std::vector<realT> single_PSD( freq.size() );
638
639 for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
640 {
641 singleLayerPSD( single_PSD, freq, m, n, i, p, fmax );
642
643 // Now add the single layer PSD to the overall PSD, weighted by Cn2
644 for( size_t j = 0; j < freq.size(); ++j )
645 {
646 PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j];
647 }
648 }
649
650 return 0;
651}
652
653template <typename realT, typename aosysT>
654template <bool parallel>
655int fourierTemporalPSD<realT, aosysT>::multiLayerPSD(
656 std::vector<realT> &PSD, std::vector<realT> &freq, realT m, realT n, int p, realT fmax )
657{
658 // PSD is zeroed every time to make sure we don't accumulate on repeated calls
659 for( size_t j = 0; j < PSD.size(); ++j )
660 PSD[j] = 0;
661
662 if( fmax == 0 )
663 {
664 fmax = 150 + 2 * fastestPeak( m, n );
665 }
666
667 return m_multiLayerPSD( PSD, freq, m, n, p, fmax, isParallel<parallel>() );
668}
669
670template <typename realT, typename aosysT>
671void fourierTemporalPSD<realT, aosysT>::makePSDGrid(
672 const std::string &dir, int mnMax, realT dFreq, realT maxFreq, realT fmax )
673{
674 std::vector<realT> freq;
675
676 std::vector<sigproc::fourierModeDef> spf;
677
678 std::string fn;
679
680 sigproc::makeFourierModeFreqs_Rect( spf, 2 * mnMax );
681
682 // Calculate number of samples, and make sure we get to at least maxFreq
683 int N = (int)( maxFreq / dFreq );
684 if( N * dFreq < maxFreq )
685 N += 1;
686
687 /*** Dump Params to file ***/
688 mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
689
690 std::ofstream fout;
691 fn = dir + '/' + "params.txt";
692 fout.open( fn );
693
694 fout << "#---------------------------\n";
695 m_aosys->dumpAOSystem( fout );
696 fout << "#---------------------------\n";
697 fout << "# PSD Grid Parameters\n";
698 fout << "# absTol " << _absTol << '\n';
699 fout << "# relTol " << _relTol << '\n';
700 fout << "# useBasis " << _useBasis << '\n';
701 fout << "# makePSDGrid call:\n";
702 fout << "# mnMax = " << mnMax << '\n';
703 fout << "# dFreq = " << dFreq << '\n';
704 fout << "# maxFreq = " << maxFreq << '\n';
705 fout << "# fmax = " << fmax << '\n';
706 fout << "#---------------------------\n";
707
708 fout.close();
709
710 // Make directory
711 std::string psddir = dir + "/psds";
712 mkdir( psddir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
713
714 // Create frequency scale.
715 math::vectorScale( freq, N, dFreq, 0 ); // dFreq); //offset from 0 by dFreq, so f=0 not included
716
717 fn = psddir + '/' + "freq.binv";
718
719 ioutils::writeBinVector( fn, freq );
720
721 size_t nLoops = 0.5 * spf.size();
722
723 ipc::ompLoopWatcher<> watcher( nLoops, std::cout );
724
725#pragma omp parallel
726 {
727 std::vector<realT> PSD;
728 PSD.resize( N );
729 std::string fname;
730
731 int m, n;
732
733#pragma omp for
734 for( size_t i = 0; i < nLoops; ++i )
735 {
736 m = spf[i * 2].m;
737 n = spf[i * 2].n;
738
739 if( fabs( (realT)m / m_aosys->D() ) >= m_aosys->spatialFilter_ku() ||
740 fabs( (realT)n / m_aosys->D() ) >= m_aosys->spatialFilter_kv() )
741 {
742 watcher.incrementAndOutputStatus();
743 continue;
744 }
745
746 multiLayerPSD<false>( PSD, freq, m, n, 1, fmax );
747
748 fname =
749 psddir + '/' + "psd_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
750
751 ioutils::writeBinVector( fname, PSD );
752
753 watcher.incrementAndOutputStatus();
754 }
755 }
756}
757
758template <typename realT, typename aosysT>
759int fourierTemporalPSD<realT, aosysT>::analyzePSDGrid( const std::string &subDir,
760 const std::string &psdDir,
761 int mnMax,
762 int mnCon,
763 realT gfixed,
764 int lpNc,
765 realT lpRegPrecision,
766 std::vector<realT> &mags,
767 int lifetimeTrials,
768 bool uncontrolledLifetimes,
769 bool writePSDs,
770 bool writeXfer )
771{
772
773 std::string dir = psdDir + "/" + subDir;
774
775 /*** Dump Params to file ***/
776 mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
777
778 std::ofstream fout;
779 std::string fn = dir + '/' + "params.txt";
780 fout.open( fn );
781
782 fout << "#---------------------------\n";
783 m_aosys->dumpAOSystem( fout );
784 fout << "#---------------------------\n";
785 fout << "# Analysis Parameters\n";
786 fout << "# mnMax = " << mnMax << '\n';
787 fout << "# mnCon = " << mnCon << '\n';
788 fout << "# lpNc = " << lpNc << '\n';
789 fout << "# mags = ";
790 for( size_t i = 0; i < mags.size() - 1; ++i )
791 fout << mags[i] << ",";
792 fout << mags[mags.size() - 1] << '\n';
793 fout << "# lifetimeTrials = " << lifetimeTrials << '\n';
794 fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n';
795 fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n';
796 fout << "# writeXfer = " << std::boolalpha << writeXfer << '\n';
797
798 fout.close();
799
800 //**** Calculating A Variance Map ****//
801
802 realT fs = 1.0 / m_aosys->tauWFS();
803 realT tauWFS = m_aosys->tauWFS();
804 realT deltaTau = m_aosys->deltaTau();
805
806 std::vector<sigproc::fourierModeDef> fms;
807
808 sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax );
809 size_t nModes = 0.5 * fms.size();
810
811 Eigen::Array<realT, -1, -1> gains, vars, speckleLifetimes, gains_lp, vars_lp, speckleLifetimes_lp;
812
813 gains.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
814 vars.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
815 speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
816
817 gains( mnMax, mnMax ) = 0;
818 vars( mnMax, mnMax ) = 0;
819 speckleLifetimes( mnMax, mnMax ) = 0;
820
821 gains_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
822 vars_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
823 speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
824
825 gains_lp( mnMax, mnMax ) = 0;
826 vars_lp( mnMax, mnMax ) = 0;
827 speckleLifetimes_lp( mnMax, mnMax ) = 0;
828
829 bool doLP = false;
830 if( lpNc > 1 )
831 doLP = true;
832 Eigen::Array<realT, -1, -1> lpC;
833
834 if( doLP )
835 {
836 lpC.resize( nModes, lpNc );
837 lpC.setZero();
838 }
839
840 std::vector<realT> S_si, S_lp;
841
842 if( writePSDs )
843 {
844 for( size_t s = 0; s < mags.size(); ++s )
845 {
846 std::string psdOutDir = dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_si";
847 mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
848
849 if( doLP )
850 {
851 psdOutDir = dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_lp";
852 mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
853 }
854 }
855 }
856
857 m_aosys->beta_p( 1, 1 );
858
859 ipc::ompLoopWatcher<> watcher( nModes * mags.size(), std::cout );
860
861 for( size_t s = 0; s < mags.size(); ++s )
862 {
863 m_aosys->starMag( mags[s] );
864
865 // In non-parallel space calculate OG=Strehl
866
867 realT opticalGain{ 1.0 };
868
869 if( m_strehlOG )
870 {
871 // Iterative optical gain
872 /// \todo need upstream NCP and CP NCP and NCP NCP
873 realT ncp = m_aosys->ncp_wfe(); // save ncp and then set it to zero for this part.
874 m_aosys->ncp_wfe( 0 );
875
876 realT lam_sci = m_aosys->lam_sci();
877 m_aosys->lam_sci( m_aosys->lam_wfs() );
878
879 realT S = m_aosys->strehl();
880 std::cerr << S << "\n";
881
882 for( int s = 0; s < 4; ++s )
883 {
884 m_aosys->opticalGain( S );
885 m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
886 S = m_aosys->strehl();
887 std::cerr << S << "\n";
888 }
889
890 opticalGain = S;
891
892 m_aosys->lam_sci( lam_sci );
893 m_aosys->ncp_wfe( ncp );
894 }
895
896 m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
897 realT strehl = m_aosys->strehl();
898
899#pragma omp parallel
900 {
901 realT localMag = mags[s];
902
903 realT var0;
904
905 realT gopt, var;
906
907 realT gopt_lp, var_lp;
908
909 std::vector<realT> tfreq; // The frequency scale of the PSDs
910 std::vector<realT> tPSDp; // The open-loop turbulence PSD for a Fourier mode
911 std::vector<realT>
912 tPSDpPOL; // The pseudo-open-loop turbulence PSD for a Fourier mode, with optical gain effects included
913
914 std::vector<realT> tfreqHF; // The above-Nyquist frequencies, saved if outputing the PSDS.
915 std::vector<realT> tPSDpHF; // The above-Nyquist component of the open-loop PSD, saved if outputing the
916 // PSDs.
917
918 //**< Get the frequency grid, and nyquist limit it to f_s/2
919 getGridPSD( tfreq, tPSDp, psdDir, 0, 1 ); // To get the freq grid
920
921 size_t imax = 0;
922 while( tfreq[imax] <= 0.5 * fs )
923 {
924 ++imax;
925 if( imax > tfreq.size() - 1 )
926 break;
927 }
928
929 if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
930 {
931 ++imax;
932 }
933
934 if( writePSDs )
935 {
936 tfreqHF.assign( tfreq.begin(), tfreq.end() );
937 }
938
939 tfreq.erase( tfreq.begin() + imax, tfreq.end() );
940
941 tPSDpPOL.resize( tfreq.size() ); // pre=allocate
942 //**>
943
944 std::vector<realT> tPSDn; // The open-loop WFS noise PSD
945 tPSDn.resize( tfreq.size() );
946
947 //**< Setup the controllers
948 mx::AO::analysis::clAOLinearPredictor<realT> tflp;
949 tflp.m_precision0 = lpRegPrecision;
950 /*tflp.m_min_sc0 = 0;
951 tflp.m_max_sc0 = 1000;*/
952
953 mx::AO::analysis::clGainOpt<realT> go_si( tauWFS, deltaTau );
954 mx::AO::analysis::clGainOpt<realT> go_lp( tauWFS, deltaTau );
955
956 go_si.f( tfreq );
957 go_lp.f( tfreq );
958
959 realT gmax = 0;
960 realT gmax_lp = 0;
961 //**>
962
963 int m, n;
964
965 //**< For use in lifetime calculations
966 sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
967 std::vector<std::complex<realT>> ETFxn;
968 std::vector<std::complex<realT>> NTFxn;
969
970 if( lifetimeTrials > 0 )
971 {
972 ETFxn.resize( tfreq.size() );
973 NTFxn.resize( tfreq.size() );
974
975 if( writeXfer )
976 {
977 std::string tfOutFile = dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_si/";
978 ioutils::createDirectories( tfOutFile );
979 }
980
981 if( doLP )
982 {
983 if( writeXfer )
984 {
985 std::string tfOutFile = dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_lp/";
986 ioutils::createDirectories( tfOutFile );
987 }
988 }
989 }
990
991//**>
992
993/*#pragma omp critical
994std::cerr << __FILE__ << " " << __LINE__ << "\n";
995*/
996// want to schedule dynamic with small chunks so maximal processor usage,
997// otherwise we can end up with a small number of cores being used at the end
998#pragma omp for schedule( dynamic, 5 )
999 for( size_t i = 0; i < nModes; ++i )
1000 {
1001 // Determine the spatial frequency at this step
1002 m = fms[2 * i].m;
1003 n = fms[2 * i].n;
1004
1005 if( fabs( (realT)m / m_aosys->D() ) >= m_aosys->spatialFilter_ku() ||
1006 fabs( (realT)n / m_aosys->D() ) >= m_aosys->spatialFilter_kv() )
1007 {
1008 gains( mnMax + m, mnMax + n ) = 0;
1009 gains( mnMax - m, mnMax - n ) = 0;
1010
1011 gains_lp( mnMax + m, mnMax + n ) = 0;
1012 gains_lp( mnMax - m, mnMax - n ) = 0;
1013
1014 vars( mnMax + m, mnMax + n ) = 0;
1015 vars( mnMax - m, mnMax - n ) = 0;
1016
1017 vars_lp( mnMax + m, mnMax + n ) = 0;
1018 vars_lp( mnMax - m, mnMax - n ) = 0;
1019 speckleLifetimes( mnMax + m, mnMax + n ) = 0;
1020 speckleLifetimes( mnMax - m, mnMax - n ) = 0;
1021 speckleLifetimes_lp( mnMax + m, mnMax + n ) = 0;
1022 speckleLifetimes_lp( mnMax - m, mnMax - n ) = 0;
1023 }
1024 else
1025 {
1026
1027 realT k = sqrt( m * m + n * n ) / m_aosys->D();
1028
1029 //**< Get the open-loop turb. PSD
1030 getGridPSD( tPSDp, psdDir, m, n );
1031
1032 // Get integral of entire open-loop PSD
1033 var0 = sigproc::psdVar( tfreq, tPSDp );
1034
1035 if( writePSDs )
1036 {
1037 tPSDpHF.assign( tPSDp.begin() + imax, tPSDp.end() );
1038 }
1039
1040 // erase points above Nyquist limit
1041 tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() );
1042
1043 // And now determine the variance which has been erased.
1044 // limVar is the out-of-band variance, which we add back in for completeness
1045 realT limVar = 0; // var0 - sigproc::psdVar( tfreq, tPSDp);
1046
1047 // And construct the POL psd
1048 if( m_uncorrectedOG )
1049 {
1050 for( size_t n = 0; n < tPSDp.size(); ++n )
1051 {
1052 tPSDpPOL[n] = tPSDp[n] * pow( opticalGain, 2 );
1053 }
1054 }
1055 else
1056 {
1057 for( size_t n = 0; n < tPSDp.size(); ++n )
1058 {
1059 tPSDpPOL[n] = tPSDp[n];
1060 }
1061 }
1062 //**>
1063
1064 //**< Determine if we're inside the hardwarecontrol limit
1065 bool inside = false;
1066
1067 if( m_aosys->circularLimit() )
1068 {
1069 if( m * m + n * n <= mnCon * mnCon )
1070 inside = true;
1071 }
1072 else
1073 {
1074 if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
1075 inside = true;
1076 }
1077 //**>
1078
1079 /* This is to select specific points for troubleshooting*/
1080 // if( !( ( m ==-2 && n == 84 ) || (m==-2 && n == 2800)) ) inside = false;
1081
1082 if( inside )
1083 {
1084 // Get the WFS noise PSD (which is already resized to match tfreq)
1085 wfsNoisePSD<realT>( tPSDn,
1086 m_aosys->beta_p( m, n ) / sqrt( opticalGain ),
1087 m_aosys->Fg( localMag ),
1088 tauWFS,
1089 m_aosys->npix_wfs( (size_t)0 ),
1090 m_aosys->Fbg( (size_t)0 ),
1091 m_aosys->ron_wfs( (size_t)0 ) );
1092
1093 gmax = 0;
1094 if( gfixed > 0 )
1095 {
1096 gopt = gfixed;
1097 var = go_si.clVariance( tPSDp, tPSDn, gopt );
1098 }
1099 else
1100 {
1101 // Calculate gain using the POL PSD
1102 gopt = go_si.optGainOpenLoop( var, tPSDpPOL, tPSDn, gmax, true );
1103
1104 if( m_uncorrectedOG )
1105 {
1106 gopt *= opticalGain;
1107 }
1108
1109 // But use the variance from the actual POL
1110 var = go_si.clVariance( tPSDp, tPSDn, gopt );
1111
1112 // Output gain curve for this mode (for troubleshooting)
1113 /*#pragma omp critical
1114 {
1115 std::string foutn = "gcurve_";
1116 foutn += std::to_string(m) + "_" + std::to_string(n) + ".dat";
1117
1118 std::ofstream foutf(foutn);
1119
1120 for(size_t n = 0; n < 10000; ++n)
1121 {
1122 realT gg = (1.0*n)/10000.;
1123 foutf << gg << " " << go_si.clVariance(tPSDp, tPSDn, gg) << "\n";
1124 }
1125
1126 foutf.close();
1127
1128 std::cerr << "\n" << gmax << " " << gopt << " " << var << " " << go_si.clVariance(tPSDp,
1129 tPSDn, 0.64) << "\n";
1130 }*/
1131 }
1132
1133 var += limVar;
1134
1135 if( doLP )
1136 {
1137 realT min_sc;
1138 int rv = tflp.regularizeCoefficients( gmax_lp,
1139 gopt_lp,
1140 var_lp,
1141 min_sc,
1142 go_lp,
1143 tPSDpPOL,
1144 tPSDn,
1145 lpNc );
1146
1147 if( rv < 0 )
1148 {
1149 std::cerr
1150 << "fourierTemporalPSD::analyzePSDGrid: regularizeCoefficients returned error ";
1151 std::cerr << rv << ' ';
1152 std::cerr << __FILE__ << ' ' << __LINE__ << '\n';
1153 }
1154
1155 for( int n = 0; n < lpNc; ++n )
1156 {
1157 lpC( i, n ) = go_lp.a( n );
1158 }
1159
1160 if( m_uncorrectedOG )
1161 {
1162 go_lp.aScale( opticalGain );
1163 go_lp.bScale( opticalGain );
1164 gopt_lp *= opticalGain;
1165 }
1166
1167 var_lp = go_lp.clVariance( tPSDp, tPSDn, gopt_lp );
1168 var_lp += limVar;
1169 }
1170 else
1171 {
1172 gopt_lp = 0;
1173 }
1174 }
1175 else
1176 {
1177 // Zero the noise PSD
1178 tPSDn.assign( tPSDn.size(), 0.0 );
1179
1180 gopt = 0;
1181 var = var0;
1182 var_lp = var0;
1183 gopt_lp = 0;
1184 go_lp.a( std::vector<realT>( { 1 } ) );
1185 go_lp.b( std::vector<realT>( { 1 } ) );
1186 }
1187
1188 //**< Determine if closed-loop is making a difference:
1189
1190 if( gopt > 0 && var > var0 )
1191 {
1192 gopt = 0;
1193 var = var0;
1194 }
1195
1196 if( gopt_lp > gopt && var_lp > var )
1197 {
1198 //Set LP to SI (or off if SI is off)
1199 gopt_lp = gopt;
1200 var_lp = var;
1201 go_lp.a( std::vector<realT>( { 1 } ) );
1202 go_lp.b( std::vector<realT>( { 1 } ) );
1203 }
1204 //**>
1205
1206 //**< Fill in the gain and variance maps
1207 gains( mnMax + m, mnMax + n ) = gopt;
1208 gains( mnMax - m, mnMax - n ) = gopt;
1209
1210 gains_lp( mnMax + m, mnMax + n ) = gopt_lp;
1211 gains_lp( mnMax - m, mnMax - n ) = gopt_lp;
1212
1213 vars( mnMax + m, mnMax + n ) = var;
1214 vars( mnMax - m, mnMax - n ) = var;
1215
1216 vars_lp( mnMax + m, mnMax + n ) = var_lp;
1217 vars_lp( mnMax - m, mnMax - n ) = var_lp;
1218 //**>
1219
1220 //**< Calculate Speckle Lifetimes
1221 if( ( lifetimeTrials > 0 || writeXfer ) && ( uncontrolledLifetimes || inside ) )
1222 {
1223 std::vector<realT> spfreq, sppsd;
1224
1225 if( gopt > 0 )
1226 {
1227 for( size_t i = 0; i < tfreq.size(); ++i )
1228 {
1229 ETFxn[i] = go_si.clETF( i, gopt );
1230 NTFxn[i] = go_si.clNTF( i, gopt );
1231 }
1232 }
1233 else
1234 {
1235 for( size_t i = 0; i < tfreq.size(); ++i )
1236 {
1237 ETFxn[i] = 1;
1238 NTFxn[i] = 0;
1239 }
1240 }
1241
1242 if( writeXfer )
1243 {
1244 std::string tfOutFile =
1245 dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_si/";
1246
1247 std::string etfOutFile = tfOutFile + "etf_" + ioutils::convertToString( m ) + '_' +
1248 ioutils::convertToString( n ) + ".binv";
1249 ioutils::writeBinVector( etfOutFile, ETFxn );
1250
1251 std::string ntfOutFile = tfOutFile + "ntf_" + ioutils::convertToString( m ) + '_' +
1252 ioutils::convertToString( n ) + ".binv";
1253 ioutils::writeBinVector( ntfOutFile, NTFxn );
1254
1255 if( i == 0 ) // Write freq on the first one
1256 {
1257 std::string fOutFile = tfOutFile + "freq.binv";
1258 ioutils::writeBinVector( fOutFile, tfreq );
1259 }
1260 }
1261
1262 if( lifetimeTrials > 0 )
1263 {
1264 speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1265 realT spvar = sigproc::psdVar( spfreq, sppsd );
1266
1267 realT splifeT = 100.0;
1268 realT error;
1269
1270 realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1271
1272 speckleLifetimes( mnMax + m, mnMax + n ) = tau;
1273 speckleLifetimes( mnMax - m, mnMax - n ) = tau;
1274 }
1275
1276 if( doLP )
1277 {
1278 if( gopt_lp > 0 )
1279 {
1280 for( size_t i = 0; i < tfreq.size(); ++i )
1281 {
1282 ETFxn[i] = go_lp.clETF( i, gopt_lp );
1283 NTFxn[i] = go_lp.clNTF( i, gopt_lp );
1284 }
1285 }
1286 else
1287 {
1288 for( size_t i = 0; i < tfreq.size(); ++i )
1289 {
1290 ETFxn[i] = 1;
1291 NTFxn[i] = 0;
1292 }
1293 }
1294
1295 if( writeXfer )
1296 {
1297 std::string tfOutFile =
1298 dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_lp/";
1299
1300 std::string etfOutFile = tfOutFile + "etf_" + ioutils::convertToString( m ) + '_' +
1301 ioutils::convertToString( n ) + ".binv";
1302 ioutils::writeBinVector( etfOutFile, ETFxn );
1303
1304 std::string ntfOutFile = tfOutFile + "ntf_" + ioutils::convertToString( m ) + '_' +
1305 ioutils::convertToString( n ) + ".binv";
1306 ioutils::writeBinVector( ntfOutFile, NTFxn );
1307
1308 if( i == 0 ) // Write freq on the first one
1309 {
1310 std::string fOutFile = tfOutFile + "freq.binv";
1311 ioutils::writeBinVector( fOutFile, tfreq );
1312 }
1313 }
1314
1315 if( lifetimeTrials > 0 )
1316 {
1317 speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1318 realT spvar = sigproc::psdVar( spfreq, sppsd );
1319
1320 realT splifeT = 100.0;
1321 realT error;
1322
1323 realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1324
1325 speckleLifetimes_lp( mnMax + m, mnMax + n ) = tau;
1326 speckleLifetimes_lp( mnMax - m, mnMax - n ) = tau;
1327 }
1328 } // if(doLP)
1329 } // if( (lifetimeTrials > 0 || writeXfer) && ( uncontrolledLifetimes || inside ))
1330 //**>
1331
1332 // Calculate the controlled PSDs and output
1333 if( writePSDs )
1334 {
1335 std::string psdOutFile =
1336 dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_si/";
1337 psdOutFile +=
1338 "psd_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1339
1340 std::vector<realT> psdOut( tPSDp.size() + tPSDpHF.size() );
1341
1342 // Calculate the output PSD if gains are applied
1343 if( gopt > 0 )
1344 {
1345 realT ETF, NTF;
1346
1347 for( size_t i = 0; i < tfreq.size(); ++i )
1348 {
1349 go_si.clTF2( ETF, NTF, i, gopt );
1350 psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
1351 }
1352
1353 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1354 {
1355 psdOut[tfreq.size() + i] = tPSDpHF[i];
1356 }
1357 }
1358 else // otherwise just copy
1359 {
1360 for( size_t i = 0; i < tfreq.size(); ++i )
1361 {
1362 psdOut[i] = tPSDp[i];
1363 }
1364
1365 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1366 {
1367 psdOut[tfreq.size() + i] = tPSDpHF[i];
1368 }
1369 }
1370
1371 ioutils::writeBinVector( psdOutFile, psdOut );
1372
1373 if( i == 0 ) // Write freq on the first one
1374 {
1375 psdOutFile =
1376 dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_si/freq.binv";
1377 ioutils::writeBinVector( psdOutFile, tfreqHF );
1378 }
1379
1380 if( doLP )
1381 {
1382 psdOutFile = dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_lp/";
1383 psdOutFile +=
1384 "psd_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1385
1386 // Calculate the output PSD if gains are applied
1387 if( gopt_lp > 0 )
1388 {
1389 realT ETF, NTF;
1390
1391 for( size_t i = 0; i < tfreq.size(); ++i )
1392 {
1393 go_lp.clTF2( ETF, NTF, i, gopt_lp );
1394 psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
1395 }
1396 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1397 {
1398 psdOut[tfreq.size() + i] = tPSDpHF[i];
1399 }
1400 }
1401 else // otherwise just copy
1402 {
1403 for( size_t i = 0; i < tfreq.size(); ++i )
1404 {
1405 psdOut[i] = tPSDp[i];
1406 }
1407
1408 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1409 {
1410 psdOut[tfreq.size() + i] = tPSDpHF[i];
1411 }
1412 }
1413
1414 ioutils::writeBinVector( psdOutFile, psdOut );
1415
1416 if( i == 0 )
1417 {
1418 psdOutFile =
1419 dir + "/" + "outputPSDs_" + ioutils::convertToString( mags[s] ) + "_lp/freq.binv";
1420 ioutils::writeBinVector( psdOutFile, tfreq );
1421 }
1422 }
1423 }
1424 }
1425 watcher.incrementAndOutputStatus();
1426
1427 } // omp for i..nModes
1428 } // omp Parallel
1429
1430 Eigen::Array<realT, -1, -1> cim;
1431
1432 fits::fitsFile<realT> ff;
1433 std::string fn = dir + "/gainmap_" + ioutils::convertToString( mags[s] ) + "_si.fits";
1434 ff.write( fn, gains );
1435
1436 fn = dir + "/varmap_" + ioutils::convertToString( mags[s] ) + "_si.fits";
1437 ff.write( fn, vars );
1438
1439 cim = vars;
1440
1441 realT Ssi = exp( -1 * cim.sum() );
1442 S_si.push_back( strehl );
1443 cim /= strehl;
1444
1445 fn = dir + "/contrast_" + ioutils::convertToString( mags[s] ) + "_si.fits";
1446 ff.write( fn, cim );
1447
1448 if( lifetimeTrials > 0 )
1449 {
1450 fn = dir + "/speckleLifetimes_" + ioutils::convertToString( mags[s] ) + "_si.fits";
1451 ff.write( fn, speckleLifetimes );
1452 }
1453
1454 if( doLP )
1455 {
1456 fn = dir + "/gainmap_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
1457 ff.write( fn, gains_lp );
1458
1459 fn = dir + "/lpcmap_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
1460 ff.write( fn, lpC );
1461
1462 fn = dir + "/varmap_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
1463 ff.write( fn, vars_lp );
1464
1465 cim = vars_lp;
1466
1467 // Scale Strehl by the ratio of total variance
1468 realT Slp = strehl * exp( -1 * cim.sum() ) /
1469 Ssi; // This is a hack until we do a real fitting error calculation or something
1470 S_lp.push_back( Slp );
1471 cim /= Slp;
1472
1473 fn = dir + "/contrast_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
1474 ff.write( fn, cim );
1475
1476 if( lifetimeTrials > 0 )
1477 {
1478 fn = dir + "/speckleLifetimes_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
1479 ff.write( fn, speckleLifetimes_lp );
1480 }
1481 }
1482
1483 } // s (mag)
1484
1485 fn = dir + "/strehl_si.txt";
1486 fout.open( fn );
1487 for( size_t i = 0; i < mags.size(); ++i )
1488 {
1489 fout << mags[i] << " " << S_si[i] << "\n";
1490 }
1491
1492 fout.close();
1493
1494 if( doLP )
1495 {
1496 fn = dir + "/strehl_lp.txt";
1497 fout.open( fn );
1498 for( size_t i = 0; i < mags.size(); ++i )
1499 {
1500 fout << mags[i] << " " << S_lp[i] << "\n";
1501 }
1502
1503 fout.close();
1504 }
1505
1506 return 0;
1507}
1508
1509template <typename realT, typename aosysT>
1510int fourierTemporalPSD<realT, aosysT>::intensityPSD(
1511 const std::string &subDir, // sub-directory of psdDir which contains the controlled system results,
1512 // and where the lifetimes will be written.
1513 const std::string &psdDir, // directory containing the PSDS
1514 const std::string &CvdPath, // path to the covariance decomposition
1515 int mnMax,
1516 int mnCon,
1517 std::vector<realT> &mags,
1518 int lifetimeTrials,
1519 bool writePSDs )
1520{
1521
1522 std::string dir = psdDir + "/" + subDir;
1523
1524 /*** Dump Params to file ***/
1525 mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
1526
1527 std::ofstream fout;
1528 std::string fn = dir + '/' + "splife_params.txt";
1529 fout.open( fn );
1530
1531 fout << "#---------------------------\n";
1532 m_aosys->dumpAOSystem( fout );
1533 fout << "#---------------------------\n";
1534 fout << "# Analysis Parameters\n";
1535 fout << "# mnMax = " << mnMax << '\n';
1536 fout << "# mnCon = " << mnCon << '\n';
1537 fout << "# mags = ";
1538 for( size_t i = 0; i < mags.size() - 1; ++i )
1539 fout << mags[i] << ",";
1540 fout << mags[mags.size() - 1] << '\n';
1541 fout << "# lifetimeTrials = " << lifetimeTrials << '\n';
1542 // fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n';
1543 fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n';
1544
1545 fout.close();
1546
1547 realT fs = 1.0 / m_aosys->tauWFS();
1548 realT tauWFS = m_aosys->tauWFS();
1549 realT deltaTau = m_aosys->deltaTau();
1550
1551 //** Get the Fourier Grid
1552 std::vector<sigproc::fourierModeDef> fms;
1553
1554 sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax );
1555 size_t nModes = 0.5 * fms.size();
1556 std::cerr << "nModes: " << nModes << " (" << fms.size() << ")\n";
1557
1558 Eigen::Array<realT, -1, -1> speckleLifetimes;
1559 Eigen::Array<realT, -1, -1> speckleLifetimes_lp;
1560
1561 speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1562 speckleLifetimes( mnMax, mnMax ) = 0;
1563
1564 speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1565 speckleLifetimes_lp( mnMax, mnMax ) = 0;
1566
1567 /*********************************************************************/
1568 // 0) Get the frequency grid, and nyquist limit it to f_s/2
1569 /*********************************************************************/
1570
1571 std::vector<realT> tfreq;
1572 std::vector<realT> tPSDp; // The open-loop OPD PSD
1573 std::vector<realT> tPSDn; // The open-loop WFS noise PSD
1574 std::vector<complexT> tETF;
1575 std::vector<complexT> tNTF;
1576
1577 if( getGridFreq( tfreq, psdDir ) < 0 )
1578 return -1;
1579
1580 size_t imax = 0;
1581 while( tfreq[imax] <= 0.5 * fs )
1582 {
1583 ++imax;
1584 if( imax > tfreq.size() - 1 )
1585 break;
1586 }
1587
1588 if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
1589 ++imax;
1590
1591 tfreq.erase( tfreq.begin() + imax, tfreq.end() );
1592
1593 // Now allocate memory
1594 tPSDn.resize( tfreq.size() );
1595 std::vector<std::vector<realT>> sqrtOPDPSD;
1596 sqrtOPDPSD.resize( nModes );
1597
1598 std::vector<std::vector<realT>> opdPSD;
1599 opdPSD.resize( nModes );
1600
1601 std::vector<realT> psd2sided;
1602
1603 // Store the mode variance for later normalization
1604 std::vector<realT> modeVar;
1605 modeVar.resize( nModes );
1606
1607 /*********************************************************************/
1608 // 1) Read in each PSD, and load it into the array in FFT order
1609 /*********************************************************************/
1610
1611 for( size_t i = 0; i < nModes; ++i )
1612 {
1613 // Determine the spatial frequency at this step
1614 int m = fms[2 * i].m;
1615 int n = fms[2 * i].n;
1616
1617 //**< Get the open-loop turb. PSD
1618 if( getGridPSD( tPSDp, psdDir, m, n ) < 0 )
1619 return -1;
1620 tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() ); // Nyquist limit
1621 modeVar[i] = sigproc::psdVar( tfreq, tPSDp );
1622
1623 // And now normalize
1624 sigproc::normPSD( tPSDp, tfreq, 1.0 ); // Normalize
1625 sigproc::augment1SidedPSD( psd2sided, tPSDp, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
1626
1627 opdPSD[i].resize( psd2sided.size() );
1628 sqrtOPDPSD[i].resize( psd2sided.size() );
1629
1630 for( size_t j = 0; j < psd2sided.size(); ++j )
1631 {
1632 opdPSD[i][j] = psd2sided[j] * modeVar[i];
1633 sqrtOPDPSD[i][j] = sqrt( psd2sided[j] ); // Store the square-root for later
1634 }
1635 }
1636
1637 size_t sz2Sided = psd2sided.size();
1638
1639 std::vector<realT> freq2sided;
1640 freq2sided.resize( sz2Sided );
1641 sigproc::augment1SidedPSDFreq( freq2sided, tfreq );
1642
1643 tPSDp.resize( tfreq.size() );
1644 tETF.resize( tfreq.size() );
1645 tNTF.resize( tfreq.size() );
1646
1647 std::vector<std::vector<realT>> sqrtNPSD;
1648 sqrtNPSD.resize( nModes );
1649
1650 std::vector<realT> noiseVar;
1651 noiseVar.resize( nModes );
1652
1653 std::vector<std::vector<complexT>> ETFsi;
1654 std::vector<std::vector<complexT>> ETFlp;
1655 ETFsi.resize( nModes );
1656 ETFlp.resize( nModes );
1657
1658 std::vector<std::vector<complexT>> NTFsi;
1659 std::vector<std::vector<complexT>> NTFlp;
1660 NTFsi.resize( nModes );
1661 NTFlp.resize( nModes );
1662
1663 std::string tfInFile;
1664 std::string etfInFile;
1665 std::string ntfInFile;
1666
1667 improc::eigenImage<realT> Cvd;
1668 fits::fitsFile<realT> ff;
1669 ff.read( Cvd, CvdPath );
1670
1671 std::vector<std::complex<realT>> tPSDc, psd2sidedc;
1672
1673 /*********************************************************************/
1674 // 2) Analyze each star magnitude
1675 /*********************************************************************/
1676 ipc::ompLoopWatcher<> watcher( lifetimeTrials * mags.size(), std::cout );
1677 for( size_t s = 0; s < mags.size(); ++s )
1678 {
1679 /*********************************************************************/
1680 // 2.0) Read in the transfer functions for each mode
1681 /*********************************************************************/
1682 for( size_t i = 0; i < nModes; ++i )
1683 {
1684 // Determine the spatial frequency at this step
1685 int m = fms[2 * i].m;
1686 int n = fms[2 * i].n;
1687
1688 //**< Determine if we're inside the hardwarecontrol limit
1689 bool inside = false;
1690
1691 if( m_aosys->circularLimit() )
1692 {
1693 if( m * m + n * n <= mnCon * mnCon )
1694 inside = true;
1695 }
1696 else
1697 {
1698 if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
1699 inside = true;
1700 }
1701
1702 // Get the WFS noise PSD (which is already resized to match tfreq)
1703 wfsNoisePSD<realT>( tPSDn,
1704 m_aosys->beta_p( m, n ),
1705 m_aosys->Fg( mags[s] ),
1706 tauWFS,
1707 m_aosys->npix_wfs( (size_t)0 ),
1708 m_aosys->Fbg( (size_t)0 ),
1709 m_aosys->ron_wfs( (size_t)0 ) );
1710 sigproc::augment1SidedPSD( psd2sided, tPSDn, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
1711
1712 // Pre-calculate the variance of the noise for later use
1713 noiseVar[i] = sigproc::psdVar( tfreq, tPSDn );
1714
1715 sqrtNPSD[i].resize( psd2sided.size() );
1716 for( size_t j = 0; j < psd2sided.size(); ++j )
1717 sqrtNPSD[i][j] = sqrt( psd2sided[j] );
1718
1719 ETFsi[i].resize( sz2Sided );
1720 ETFlp[i].resize( sz2Sided );
1721 NTFsi[i].resize( sz2Sided );
1722 NTFlp[i].resize( sz2Sided );
1723
1724 if( inside )
1725 {
1726 tfInFile = dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_si/";
1727
1728 etfInFile =
1729 tfInFile + "etf_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1730 ioutils::readBinVector( tPSDc, etfInFile );
1731 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1732 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1733 ETFsi[i][j] = psd2sidedc[j];
1734
1735 ntfInFile =
1736 tfInFile + "ntf_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1737 ioutils::readBinVector( tPSDc, ntfInFile );
1738 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1739 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1740 NTFsi[i][j] = psd2sidedc[j];
1741
1742 tfInFile = dir + "/" + "outputTF_" + ioutils::convertToString( mags[s] ) + "_lp/";
1743
1744 etfInFile =
1745 tfInFile + "etf_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1746 ioutils::readBinVector( tPSDc, etfInFile );
1747 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1748 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1749 ETFlp[i][j] = psd2sidedc[j];
1750
1751 ntfInFile =
1752 tfInFile + "ntf_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
1753 ioutils::readBinVector( tPSDc, ntfInFile );
1754 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1755 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1756 NTFlp[i][j] = psd2sidedc[j];
1757 }
1758 else
1759 {
1760 for( int q = 0; q < ETFsi.size(); ++q )
1761 {
1762 ETFsi[i][q] = 1;
1763 NTFsi[i][q] = 0;
1764 ETFlp[i][q] = 1;
1765 NTFlp[i][q] = 0;
1766 }
1767 }
1768 }
1769
1770 sigproc::averagePeriodogram<realT> tavgPgram(
1771 sz2Sided / 1.,
1772 1 / fs ); // this is just to get the size right, per-thread instances below
1773 std::vector<std::vector<realT>> spPSDs;
1774 spPSDs.resize( nModes );
1775 for( size_t pp = 0; pp < spPSDs.size(); ++pp )
1776 {
1777 spPSDs[pp].resize( tavgPgram.size() );
1778 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
1779 spPSDs[pp][nn] = 0;
1780 }
1781
1782 std::vector<std::vector<realT>> spPSDslp;
1783 spPSDslp.resize( nModes );
1784 for( size_t pp = 0; pp < spPSDslp.size(); ++pp )
1785 {
1786 spPSDslp[pp].resize( tavgPgram.size() );
1787 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
1788 spPSDslp[pp][nn] = 0;
1789 }
1790
1791#pragma omp parallel
1792 {
1793 // Normally distributed random numbers
1794 math::normDistT<realT> normVar;
1795 normVar.seed();
1796
1797 // FFTs for going to Fourier domain and back to time domain.
1798 math::ft::fftT<realT, std::complex<realT>, 1, 0> fftF( sqrtOPDPSD[0].size() );
1799 math::ft::fftT<std::complex<realT>, realT, 1, 0> fftB( sqrtOPDPSD[0].size(), math::ft::dir::backward );
1800
1801 // Working memory
1802 std::vector<std::complex<realT>> tform1( sqrtOPDPSD[0].size() );
1803 std::vector<std::complex<realT>> tform2( sqrtOPDPSD[0].size() );
1804 std::vector<std::complex<realT>> Ntform1( sqrtOPDPSD[0].size() );
1805 std::vector<std::complex<realT>> Ntform2( sqrtOPDPSD[0].size() );
1806
1807 std::vector<std::complex<realT>> tform1lp( sqrtOPDPSD[0].size() );
1808 std::vector<std::complex<realT>> tform2lp( sqrtOPDPSD[0].size() );
1809 std::vector<std::complex<realT>> Ntform1lp( sqrtOPDPSD[0].size() );
1810 std::vector<std::complex<realT>> Ntform2lp( sqrtOPDPSD[0].size() );
1811
1812 // OPD-PSD filter
1813 sigproc::psdFilter<realT, 1> pfilt;
1814 pfilt.psdSqrt( &sqrtOPDPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
1815
1816 // Noise-PSD filter
1817 sigproc::psdFilter<realT, 1> nfilt;
1818 nfilt.psdSqrt( &sqrtNPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
1819
1820 // The h time-series
1821 std::vector<std::vector<realT>> hts;
1822 hts.resize( 2 * nModes );
1823
1824 // The correlated h time-series
1825 std::vector<std::vector<realT>> htsCorr;
1826 htsCorr.resize( 2 * nModes );
1827
1828 for( size_t pp = 0; pp < hts.size(); ++pp )
1829 {
1830 hts[pp].resize( sqrtOPDPSD[0].size() );
1831 htsCorr[pp].resize( sqrtOPDPSD[0].size() );
1832 }
1833
1834 // The noise time-serieses
1835 std::vector<realT> N_n;
1836 N_n.resize( sz2Sided );
1837
1838 std::vector<realT> N_nm;
1839 N_nm.resize( sz2Sided );
1840
1841 // Periodogram averager
1842 sigproc::averagePeriodogram<realT> avgPgram( sz2Sided / 1., 1 / fs ); //, 0, 1);
1843 avgPgram.win( sigproc::window::hann );
1844
1845 // The periodogram output
1846 std::vector<realT> tpgram( avgPgram.size() );
1847
1848 // Holds the speckle time-series
1849 improc::eigenCube<realT> spTS;
1850 spTS.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
1851
1852 improc::eigenCube<realT> spTSlp;
1853 spTSlp.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
1854
1855// Here's where the big loop of n-trials should start
1856#pragma omp for
1857 for( int zz = 0; zz < lifetimeTrials; ++zz )
1858 {
1859 std::complex<realT> scale = exp( std::complex<realT>( 0, math::half_pi<realT>() ) ) /
1860 std::complex<realT>( ( tform1.size() ), 0 );
1861
1862 /*********************************************************************/
1863 // 2.1) Generate filtered noise for each mode, with temporal phase shifts at each spatial frequency
1864 /*********************************************************************/
1865 for( size_t pp = 0; pp < nModes; ++pp )
1866 {
1867 // Fill in standard normal noise
1868 for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
1869 {
1870 hts[2 * pp][nn] = normVar;
1871 }
1872
1873 // Set sqrt(PSD), just a pointer switch
1874 pfilt.psdSqrt( &sqrtOPDPSD[pp], tfreq[1] - tfreq[0] );
1875
1876 // And now filter the noise to a time-series of h
1877 pfilt( hts[2 * pp] );
1878
1879 /**/
1880 // Then construct 2nd mode with temporal shift
1881 fftF( tform1.data(), hts[2 * pp].data() );
1882
1883 // Apply the phase shift to form the 2nd time series
1884 for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
1885 tform1[nn] = tform1[nn] * scale;
1886
1887 fftB( hts[2 * pp + 1].data(), tform1.data() );
1888 /**/
1889 }
1890 //** At this point we have correlated time-series, with the correct temporal PSD, but not yet spatially
1891 // correlated
1892
1893 /*********************************************************************/
1894 // 2.2) Correlate the time-series for each mode
1895 /*********************************************************************/
1896 // #pragma omp parallel for
1897 for( size_t pp = 0; pp < hts.size(); ++pp )
1898 {
1899 for( size_t nn = 0; nn < hts[0].size(); ++nn )
1900 {
1901 htsCorr[pp][nn] = 0;
1902 }
1903
1904 for( size_t qq = 0; qq <= pp; ++qq )
1905 {
1906 realT cvd = Cvd( qq, pp );
1907 realT *d1 = htsCorr[pp].data();
1908 realT *d2 = hts[qq].data();
1909 for( size_t nn = 0; nn < hts[0].size(); ++nn )
1910 {
1911 d1[nn] += d2[nn] * cvd;
1912 }
1913 }
1914 }
1915
1916 /*
1917 for(size_t pp=0; pp < hts.size(); ++pp)
1918 {
1919 for(size_t nn=0; nn< hts[0].size(); ++nn)
1920 {
1921 htsCorr[pp][nn] = hts[pp][nn];
1922 }
1923 }*/
1924
1925 /*********************************************************************/
1926 // 2.2.a) Re-normalize b/c the correlation step above does not result in correct variances
1927 ///\todo should be able to scale the covar by r0, and possibly D
1928 /*********************************************************************/
1929 for( size_t pp = 0; pp < nModes; ++pp )
1930 {
1931 math::vectorMeanSub( htsCorr[2 * pp] );
1932 math::vectorMeanSub( htsCorr[2 * pp + 1] );
1933
1934 realT var = math::vectorVariance( htsCorr[2 * pp] );
1935 realT norm = sqrt( modeVar[pp] / var );
1936 for( size_t nn = 0; nn < htsCorr[2 * pp].size(); ++nn )
1937 htsCorr[2 * pp][nn] *= norm;
1938
1939 var = math::vectorVariance( htsCorr[2 * pp + 1] );
1940 norm = sqrt( modeVar[pp] / var );
1941 for( size_t nn = 0; nn < htsCorr[2 * pp + 1].size(); ++nn )
1942 htsCorr[2 * pp + 1][nn] *= norm;
1943 }
1944
1945 scale = std::complex<realT>( tform1.size(), 0 );
1946
1947 /*********************************************************************/
1948 // 2.3) Generate speckle intensity time-series
1949 /*********************************************************************/
1950 for( size_t pp = 0; pp < nModes; ++pp )
1951 {
1952 // Now we take them back to the FD and apply the xfer
1953 // and add the noise
1954
1955 fftF( tform1.data(), htsCorr[2 * pp].data() );
1956 fftF( tform2.data(), htsCorr[2 * pp + 1].data() );
1957
1958 // Make some noise
1959 for( int nn = 0; nn < sz2Sided; ++nn )
1960 {
1961 N_n[nn] = normVar;
1962 N_nm[nn] = normVar;
1963 }
1964
1965 // Filter it
1966 // Set sqrt(PSD), just a pointer switch
1967 pfilt.psdSqrt( &sqrtNPSD[pp], tfreq[1] - tfreq[0] );
1968 nfilt.filter( N_n );
1969 nfilt.filter( N_nm );
1970
1971 // Normalize it
1972 realT Nactvar = 0.5 * ( math::vectorVariance( N_n ) + math::vectorVariance( N_nm ) );
1973 realT norm = sqrt( noiseVar[pp] / Nactvar );
1974 for( size_t q = 0; q < N_n.size(); ++q )
1975 N_n[q] *= norm;
1976 for( size_t q = 0; q < N_nm.size(); ++q )
1977 N_nm[q] *= norm;
1978
1979 // And move them to the Fourier domain
1980 fftF( Ntform1.data(), N_n.data() );
1981 fftF( Ntform2.data(), N_nm.data() );
1982
1983 // Apply the closed loop transfers
1984 for( size_t mm = 0; mm < tform1.size(); ++mm )
1985 {
1986 // Apply the augmented ETF to two time-series
1987 tform1lp[mm] = tform1[mm] * ETFlp[pp][mm] / scale;
1988 tform2lp[mm] = tform2[mm] * ETFlp[pp][mm] / scale;
1989
1990 Ntform1lp[mm] = Ntform1[mm] * NTFlp[pp][mm] / scale;
1991 Ntform2lp[mm] = Ntform2[mm] * NTFlp[pp][mm] / scale;
1992
1993 tform1[mm] *= ETFsi[pp][mm] / scale;
1994 tform2[mm] *= ETFsi[pp][mm] / scale;
1995
1996 Ntform1[mm] *= NTFsi[pp][mm] / scale;
1997 Ntform2[mm] *= NTFsi[pp][mm] / scale;
1998 }
1999
2000 // And make the speckle TS
2001 int m = fms[2 * pp].m;
2002 int n = fms[2 * pp].n;
2003
2004 //<<<<<<<<****** Transform back to the time domain.
2005 fftB( htsCorr[2 * pp].data(), tform1.data() );
2006 fftB( htsCorr[2 * pp + 1].data(), tform2.data() );
2007 fftB( N_n.data(), Ntform1.data() );
2008 fftB( N_nm.data(), Ntform2.data() );
2009
2010 for( int i = 0; i < hts[2 * pp].size(); ++i )
2011 {
2012 realT h1 = htsCorr[2 * pp][i] + N_n[i];
2013 realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2014
2015 spTS.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2016 spTS.image( i )( mnMax - m, mnMax - n ) = spTS.image( i )( mnMax + m, mnMax + n );
2017 }
2018
2019 fftB( htsCorr[2 * pp].data(), tform1lp.data() );
2020 fftB( htsCorr[2 * pp + 1].data(), tform2lp.data() );
2021 fftB( N_n.data(), Ntform1lp.data() );
2022 fftB( N_nm.data(), Ntform2lp.data() );
2023
2024 for( int i = 0; i < hts[2 * pp].size(); ++i )
2025 {
2026 realT h1 = htsCorr[2 * pp][i] + N_n[i];
2027 realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2028
2029 spTSlp.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2030 spTSlp.image( i )( mnMax - m, mnMax - n ) = spTSlp.image( i )( mnMax + m, mnMax + n );
2031 }
2032 }
2033
2034 /*********************************************************************/
2035 // 2.5) Calculate speckle PSD for each mode
2036 /*********************************************************************/
2037 std::vector<realT> speckAmp( spTS.planes() );
2038 std::vector<realT> speckAmplp( spTS.planes() );
2039
2040 for( size_t pp = 0; pp < nModes; ++pp )
2041 {
2042 int m = fms[2 * pp].m;
2043 int n = fms[2 * pp].n;
2044
2045 realT mn = 0;
2046 realT mnlp = 0;
2047 for( int i = 0; i < spTS.planes(); ++i )
2048 {
2049 speckAmp[i] = spTS.image( i )( mnMax + m, mnMax + n );
2050 speckAmplp[i] = spTSlp.image( i )( mnMax + m, mnMax + n );
2051
2052 mn += speckAmp[i];
2053 mnlp += speckAmplp[i];
2054 }
2055 mn /= speckAmp.size();
2056 mnlp /= speckAmplp.size();
2057
2058 // mean subtract
2059 for( int i = 0; i < speckAmp.size(); ++i )
2060 speckAmp[i] -= mn;
2061 for( int i = 0; i < speckAmplp.size(); ++i )
2062 speckAmplp[i] -= mnlp;
2063
2064 // Calculate PSD of the speckle amplitude
2065 avgPgram( tpgram, speckAmp );
2066 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2067 spPSDs[pp][nn] += tpgram[nn];
2068
2069 avgPgram( tpgram, speckAmplp );
2070 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2071 spPSDslp[pp][nn] += tpgram[nn];
2072 }
2073
2074 watcher.incrementAndOutputStatus();
2075 } // for(int zz=0; zz<lifetimeTrials; ++zz)
2076 } // #pragma omp parallel
2077
2078 std::vector<realT> spFreq( spPSDs[0].size() );
2079 for( size_t nn = 0; nn < spFreq.size(); ++nn )
2080 spFreq[nn] = tavgPgram[nn];
2081
2082 improc::eigenImage<realT> taus, tauslp;
2083 taus.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2084 tauslp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2085
2086 improc::eigenCube<realT> imc, imclp;
2087 std::vector<realT> clPSD;
2088
2089 if( writePSDs )
2090 {
2091 imc.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2092 imclp.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2093 clPSD.resize( sz2Sided );
2094 }
2095
2096 sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
2097 /*********************************************************************/
2098 // 3.0) Calculate lifetimes from the PSDs
2099 /*********************************************************************/
2100 for( size_t pp = 0; pp < nModes; ++pp )
2101 {
2102 spPSDs[pp][0] = spPSDs[pp][1]; // deal with under-estimated mean.
2103 spPSDslp[pp][0] = spPSDslp[pp][1]; // deal with under-estimated mean.
2104
2105 int m = fms[2 * pp].m;
2106 int n = fms[2 * pp].n;
2107
2108 realT var;
2109 if( writePSDs ) // Have to normalize the intensity for some reason if we want to use the PSDs
2110 {
2111 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2112 {
2113 spPSDs[pp][nn] /= lifetimeTrials;
2114 }
2115
2116 for( size_t nn = 0; nn < sz2Sided; ++nn )
2117 {
2118 clPSD[nn] =
2119 opdPSD[pp][nn] * norm( ETFsi[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFsi[pp][nn] );
2120 }
2121
2122 var = sigproc::psdVar( freq2sided, clPSD );
2123
2124 realT pvar = sigproc::psdVar( spFreq, spPSDs[pp] );
2125
2126 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2127 {
2128 spPSDs[pp][nn] *= var / pvar;
2129 imc.image( nn )( mnMax + m, mnMax + n ) = spPSDs[pp][nn];
2130 imc.image( nn )( mnMax - m, mnMax - n ) = spPSDs[pp][nn];
2131 }
2132
2133 // lp
2134 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2135 {
2136 spPSDslp[pp][nn] /= lifetimeTrials;
2137 }
2138
2139 for( size_t nn = 0; nn < sz2Sided; ++nn )
2140 {
2141 clPSD[nn] =
2142 opdPSD[pp][nn] * norm( ETFlp[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFlp[pp][nn] );
2143 }
2144
2145 var = sigproc::psdVar( freq2sided, clPSD );
2146
2147 pvar = sigproc::psdVar( spFreq, spPSDslp[pp] );
2148
2149 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2150 {
2151 spPSDslp[pp][nn] *= var / pvar;
2152 imclp.image( nn )( mnMax + m, mnMax + n ) = spPSDslp[pp][nn];
2153 imclp.image( nn )( mnMax - m, mnMax - n ) = spPSDslp[pp][nn];
2154 }
2155 }
2156
2157 var = sigproc::psdVar( spFreq, spPSDs[pp] );
2158
2159 realT T = ( 1.0 / ( spFreq[1] - spFreq[0] ) ) * 10;
2160 realT error;
2161 realT tau = pvm( error, spFreq, spPSDs[pp], T ) * ( T ) / var;
2162 taus( mnMax + m, mnMax + n ) = tau;
2163 taus( mnMax - m, mnMax - n ) = tau;
2164
2165 var = sigproc::psdVar( spFreq, spPSDslp[pp] );
2166
2167 tau = pvm( error, spFreq, spPSDslp[pp], T ) * ( T ) / var;
2168 tauslp( mnMax + m, mnMax + n ) = tau;
2169 tauslp( mnMax - m, mnMax - n ) = tau;
2170 }
2171
2172 /*********************************************************************/
2173 // 4.0) Write the results to disk
2174 /*********************************************************************/
2175 fn = dir + "/speckleLifetimes_" + ioutils::convertToString( mags[s] ) + "_si.fits";
2176 ff.write( fn, taus );
2177
2178 fn = dir + "/speckleLifetimes_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
2179 ff.write( fn, tauslp );
2180
2181 if( writePSDs )
2182 {
2183 fn = dir + "/specklePSDs_" + ioutils::convertToString( mags[s] ) + "_si.fits";
2184 ff.write( fn, imc );
2185
2186 fn = dir + "/specklePSDs_" + ioutils::convertToString( mags[s] ) + "_lp.fits";
2187 ff.write( fn, imclp );
2188 }
2189
2190 } // mags
2191
2192 return 0;
2193}
2194
2195template <typename realT, typename aosysT>
2196int fourierTemporalPSD<realT, aosysT>::getGridFreq( std::vector<realT> &freq, const std::string &dir )
2197{
2198 std::string fn;
2199 fn = dir + "/psds/freq.binv";
2200 return ioutils::readBinVector( freq, fn );
2201}
2202
2203template <typename realT, typename aosysT>
2204int fourierTemporalPSD<realT, aosysT>::getGridPSD( std::vector<realT> &psd, const std::string &dir, int m, int n )
2205{
2206 std::string fn;
2207 fn = dir + "/psds/psd_" + ioutils::convertToString( m ) + '_' + ioutils::convertToString( n ) + ".binv";
2208 return ioutils::readBinVector( psd, fn );
2209}
2210
2211template <typename realT, typename aosysT>
2212int fourierTemporalPSD<realT, aosysT>::getGridPSD(
2213 std::vector<realT> &freq, std::vector<realT> &psd, const std::string &dir, int m, int n )
2214{
2215 int rv = getGridFreq( freq, dir );
2216 if( rv < 0 )
2217 return rv;
2218
2219 return getGridPSD( psd, dir, m, n );
2220}
2221
2222/// Worker function for GSL Integration for the basic sin/cos Fourier modes.
2223/** \ingroup mxAOAnalytic
2224 */
2225template <typename realT, typename aosysT>
2226realT F_basic( realT kv, void *params )
2227{
2228 fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2229
2230 realT f = Fp->m_f;
2231 realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2232
2233 realT D = Fp->m_aosys->D();
2234 realT m = Fp->m_m;
2235 realT n = Fp->m_n;
2236 int p = Fp->m_p;
2237
2238 realT ku = f / v_wind;
2239
2240 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2241 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2242
2243 realT Q1 = math::func::jinc( math::pi<realT>() * D * kp );
2244
2245 realT Q2 = math::func::jinc( math::pi<realT>() * D * kpp );
2246
2247 realT Q = ( Q1 + p * Q2 );
2248
2249 realT P =
2250 Fp->m_aosys->psd( Fp->m_aosys->atm, Fp->_layer_i, sqrt( pow( ku, 2 ) + pow( kv, 2 ) ), Fp->m_aosys->secZeta() );
2251
2252 return P * Q * Q;
2253}
2254
2255template <typename realT>
2256void turbBoilCubic(
2257 realT &a, realT &b, realT &c, realT &d, const realT &kv, const realT &f, const realT &Vu, const realT &f0, int pm )
2258{
2259 a = Vu * Vu * Vu;
2260 b = -( 3 * Vu * Vu * f + pm * f0 * f0 * f0 );
2261 c = 3 * f * f * Vu;
2262 d = -( f * f * f + pm * f0 * f0 * f0 * kv * kv );
2263}
2264
2265/// Worker function for GSL Integration for the modified Fourier modes.
2266/** \ingroup mxAOAnalytic
2267 */
2268template <typename realT, typename aosysT>
2269realT F_mod( realT kv, void *params )
2270{
2271 fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2272
2273 realT f = Fp->m_f;
2274 realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2275
2276 realT D = Fp->m_aosys->D();
2277 realT m = Fp->m_m;
2278 realT n = Fp->m_n;
2279
2280 realT f0 = Fp->m_f0;
2281
2282 realT ku;
2283 if( f0 == 0 )
2284 {
2285 ku = f / v_wind;
2286
2287 if( Fp->m_spatialFilter )
2288 {
2289 // de-rotate the spatial frequency vector back to pupil coordinates
2290 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2291 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2292 // Return if spatially filtered
2293 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2294 return 0;
2295
2296 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2297 return 0;
2298 }
2299
2300 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2301 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2302
2303 realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
2304
2305 realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2306
2307 realT QQ = 2 * ( Jp * Jp + Jm * Jm );
2308
2309 realT P = Fp->m_aosys->psd( Fp->m_aosys->atm,
2310 Fp->_layer_i,
2311 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2312 Fp->m_aosys->lam_sci(),
2313 Fp->m_aosys->lam_wfs(),
2314 Fp->m_aosys->secZeta() );
2315
2316 return P * QQ;
2317 }
2318 else
2319 {
2320 realT a, b, c, d, p, q;
2321
2322 turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, 1 );
2323 mx::math::cubicDepressed( p, q, a, b, c, d );
2324 realT t = mx::math::cubicRealRoot( p, q );
2325
2326 ku = t - b / ( 3 * a );
2327
2328 if( Fp->m_spatialFilter )
2329 {
2330 // de-rotate the spatial frequency vector back to pupil coordinates
2331 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2332 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2333 // Return if spatially filtered
2334 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2335 return 0;
2336
2337 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2338 return 0;
2339 }
2340
2341 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2342 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2343
2344 realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
2345
2346 realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2347
2348 realT QQ = 2 * ( Jp * Jp + Jm * Jm );
2349
2350 realT P1 = Fp->m_aosys->psd( Fp->m_aosys->atm,
2351 Fp->_layer_i,
2352 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2353 Fp->m_aosys->lam_sci(),
2354 Fp->m_aosys->lam_wfs(),
2355 Fp->m_aosys->secZeta() );
2356
2357 P1 *= QQ;
2358
2359 turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, -1 );
2360 mx::math::cubicDepressed( p, q, a, b, c, d );
2361 t = mx::math::cubicRealRoot( p, q );
2362
2363 ku = t - b / ( 3 * a );
2364
2365 if( Fp->m_spatialFilter )
2366 {
2367 // de-rotate the spatial frequency vector back to pupil coordinates
2368 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2369 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2370 // Return if spatially filtered
2371 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2372 return 0;
2373
2374 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2375 return 0;
2376 }
2377
2378 kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2379 kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2380
2381 Jp = math::func::jinc( math::pi<realT>() * D * kp );
2382
2383 Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2384
2385 QQ = 2 * ( Jp * Jp + Jm * Jm );
2386
2387 realT P2 = Fp->m_aosys->psd( Fp->m_aosys->atm,
2388 Fp->_layer_i,
2389 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2390 Fp->m_aosys->lam_sci(),
2391 Fp->m_aosys->lam_wfs(),
2392 Fp->m_aosys->secZeta() );
2393
2394 P2 *= QQ;
2395
2396 return 0.5*(P1 + P2);
2397 }
2398}
2399
2400
2401/*extern template
2402struct fourierTemporalPSD<float, aoSystem<float, vonKarmanSpectrum<float>, std::ostream>>;*/
2403
2404extern template struct fourierTemporalPSD<double, aoSystem<double, vonKarmanSpectrum<double>, std::ostream>>;
2405
2406/*
2407extern template
2408struct fourierTemporalPSD<long double, aoSystem<long double, vonKarmanSpectrum<long double>, std::ostream>>;
2409
2410#ifdef HASQUAD
2411extern template
2412struct fourierTemporalPSD<__float128, aoSystem<__float128, vonKarmanSpectrum<__float128>, std::ostream>>;
2413#endif
2414*/
2415
2416} // namespace analysis
2417} // namespace AO
2418} // namespace mx
2419
2420#endif // fourierTemporalPSD_hpp
aoConstants.hpp
Calculate and provide constants related to adaptive optics.
aoPSDs.hpp
Spatial power spectra used in adaptive optics.
aoSystem.hpp
Declares and defines an analytical AO system.
clAOLinearPredictor.hpp
Provides a class to manage closed loop gain linear predictor determination.
clGainOpt.hpp
Provides a class to manage closed loop gain optimization.
mx::fits::fitsFile
Class to manage interactions with a FITS file.
Definition fitsFile.hpp:52
mx::fits::fitsFile::write
int write(const dataT *im, int d1, int d2, int d3, fitsHeader *head)
Write the contents of a raw array to the FITS file.
Definition fitsFile.hpp:1329
mx::fits::fitsFile::read
int read(dataT *data)
Read the contents of the FITS file into an array.
Definition fitsFile.hpp:829
mx::improc::eigenCube
An image cube with an Eigen-like API.
Definition eigenCube.hpp:30
mx::improc::eigenCube::image
Eigen::Map< Eigen::Array< dataT, Eigen::Dynamic, Eigen::Dynamic > > image(Index n)
Returns a 2D Eigen::Eigen::Map pointed at the specified image.
Definition eigenCube.hpp:382
mx::ipc::ompLoopWatcher
A class to track the number of iterations in an OMP parallelized loop.
Definition ompLoopWatcher.hpp:86
mx::ipc::ompLoopWatcher::incrementAndOutputStatus
void incrementAndOutputStatus()
Increment and output status.
Definition ompLoopWatcher.hpp:274
mx::math::randomT::seed
void seed(typename ranengT::result_type seedval)
Set the seed of the random engine.
Definition randomT.hpp:101
mx::sigproc::averagePeriodogram
Calculate the average periodogram of a time-series.
Definition averagePeriodogram.hpp:86
mx::sigproc::averagePeriodogram::resize
int resize(size_t avgLen)
Resize the periodogram working memory, setting up the 1/2-overlapped optimum case.
Definition averagePeriodogram.hpp:331
mx::sigproc::averagePeriodogram::size
size_t size()
Return the size of the periodogram.
Definition averagePeriodogram.hpp:469
mx::sigproc::averagePeriodogram::win
std::vector< realT > & win()
Get a reference to the window vector.
Definition averagePeriodogram.hpp:375
WSZ
#define WSZ
The size of the gsl_integration workspaces.
Definition fourierCovariance.hpp:44
mx::AO::analysis::modified
@ modified
The modified Fourier basis from .
Definition fourierTemporalPSD.hpp:84
mx::AO::analysis::basic
@ basic
The basic sine and cosine Fourier modes.
Definition fourierTemporalPSD.hpp:83
mx::ioutils::writeBinVector
int writeBinVector(const std::string &fname, std::vector< dataT > &vec)
Write a BinVector file to disk.
Definition binVector.hpp:336
mx::ioutils::readBinVector
int readBinVector(std::vector< dataT > &vec, const std::string &fname)
Read a BinVector file from disk.
Definition binVector.hpp:241
mx::improc::eigenImage
Eigen::Array< scalarT, -1, -1 > eigenImage
Definition of the eigenImage type, which is an alias for Eigen::Array.
Definition eigenImage.hpp:44
mx::error_t::error
@ error
A general error has occurred.
mxError
#define mxError(esrc, ecode, expl)
This reports an mxlib specific error.
Definition mxErrorOld.hpp:72
mx::ioutils::createDirectories
int createDirectories(const std::string &path)
Create a directory or directories.
Definition fileUtils.cpp:57
mx::sigproc::fourierModeCoordinates
int fourierModeCoordinates(int &m, int &n, int &p, int i)
Calculate the (m,n,p) coordinates of a Fourier mode given its index.
Definition fourierModes.hpp:541
mx::sigproc::makeFourierModeFreqs_Rect
int makeFourierModeFreqs_Rect(std::vector< fourierModeDef > &spf, int N)
Generate a rectangular spatial frequency grid.
Definition fourierModes.hpp:592
mx::math::ft::dir::backward
@ backward
Specifies the backward transform.
mx::math::func::jinc
T jinc(const T &x)
The Jinc function.
Definition jinc.hpp:61
mx::math::six_fifths
constexpr floatT six_fifths()
Return 6/5 in the specified precision.
Definition constants.hpp:404
mx::AO::analysis::F_basic
realT F_basic(realT kv, void *params)
Worker function for GSL Integration for the basic sin/cos Fourier modes.
Definition fourierTemporalPSD.hpp:2226
mx::AO::analysis::F_mod
realT F_mod(realT kv, void *params)
Worker function for GSL Integration for the modified Fourier modes.
Definition fourierTemporalPSD.hpp:2269
MXE_INVALIDARG
#define MXE_INVALIDARG
An argument was invalid.
Definition mxErrorOld.hpp:112
mx::sigproc::augment1SidedPSD
void augment1SidedPSD(vectorTout &psdTwoSided, vectorTin &psdOneSided, bool addZeroFreq=false, typename vectorTin::value_type scale=0.5)
Augment a 1-sided PSD to standard 2-sided FFT form.
Definition psdUtils.hpp:817
mx::sigproc::augment1SidedPSDFreq
void augment1SidedPSDFreq(std::vector< T > &freqTwoSided, std::vector< T > &freqOneSided)
Augment a 1-sided frequency scale to standard FFT form.
Definition psdUtils.hpp:876
mx::sigproc::psdVar
realT psdVar(const std::vector< realT > &f, const std::vector< realT > &PSD, realT half=0.5)
Calculate the variance of a 1-D PSD.
Definition psdUtils.hpp:136
mx::sigproc::normPSD
int normPSD(std::vector< floatT > &psd, std::vector< floatT > &f, floatParamT normT, floatT fmin=std::numeric_limits< floatT >::min(), floatT fmax=std::numeric_limits< floatT >::max())
Normalize a 1-D PSD to have a given variance.
Definition psdUtils.hpp:447
mx::sigproc::window::hann
void hann(realT *filt, int N)
The Hann Window.
Definition signalWindows.hpp:179
gmax
#define gmax(A, B)
max(A,B) - larger (most +ve) of two numbers (generic) (defined in the SOFA library sofam....
Definition sofa_constants.h:106
mx::ioutils::convertToString
std::string convertToString(const typeT &value, int precision=0)
Convert a numerical value to a string.
Definition stringUtils.hpp:83
mx::math::vectorMeanSub
void vectorMeanSub(valueT *vec, size_t sz)
Subtract the mean from a vector.
Definition vectorUtils.hpp:530
mx::math::vectorScale
void vectorScale(vectorT &vec, size_t N=0, typename vectorT::value_type scale=0, typename vectorT::value_type offset=0)
Fill in a vector with a regularly spaced scale.
Definition vectorUtils.hpp:61
mx::math::vectorVariance
valueT vectorVariance(const valueT *vec, size_t sz, valueT mean)
Calculate the variance of a vector relative to a supplied mean value.
Definition vectorUtils.hpp:293
mx
The mxlib c++ namespace.
Definition mxError.hpp:40
speckleAmpPSD.hpp
Calculates the PSD of speckle intensity given a modified Fourier mode amplitude PSD.
mx::AO::analysis::speckleAmpPSD
int speckleAmpPSD(std::vector< realT > &spFreq, std::vector< realT > &spPSD, const std::vector< realT > &freq, const std::vector< realT > &fmPSD, const std::vector< std::complex< realT > > &fmXferFxn, const std::vector< realT > &nPSD, const std::vector< std::complex< realT > > &nXferFxn, int N, std::vector< realT > *vars=nullptr, std::vector< realT > *bins=nullptr, bool noPSD=false)
Calculate the PSD of the speckle intensity given the PSD of Fourier mode amplitude.
Definition speckleAmpPSD.hpp:50
mx::AO::analysis::clAOLinearPredictor
Class to manage the calculation of linear predictor coefficients for a closed-loop AO system.
Definition clAOLinearPredictor.hpp:42
mx::AO::analysis::clAOLinearPredictor::regularizeCoefficients
int regularizeCoefficients(realT &gmax_lp, realT &gopt_lp, realT &var_lp, realT &min_sc, clGainOpt< realT > &go_lp, std::vector< realT > &PSDt, std::vector< realT > &PSDn, int Nc)
Regularize the PSD and calculate the associated LP coefficients.
Definition clAOLinearPredictor.hpp:227
mx::AO::analysis::clGainOpt
A class to manage optimizing closed-loop gains.
Definition clGainOpt.hpp:64
mx::AO::analysis::clGainOpt::a
void a(const std::vector< realT > &newA)
Set the vector of IIR coefficients.
Definition clGainOpt.hpp:571
mx::AO::analysis::clGainOpt::clNTF
complexT clNTF(int fi, realT g)
Return the closed loop noise transfer function (NTF) at frequency f for gain g.
Definition clGainOpt.hpp:886
mx::AO::analysis::clGainOpt::b
void b(const std::vector< realT > &newB)
Set the vector of FIR coefficients.
Definition clGainOpt.hpp:530
mx::AO::analysis::clGainOpt::optGainOpenLoop
realT optGainOpenLoop(realT &var, const std::vector< realT > &PSDerr, const std::vector< realT > &PSDnoise, bool gridSearch)
Return the optimum closed loop gain given an open loop PSD.
Definition clGainOpt.hpp:1136
mx::AO::analysis::clGainOpt::clVariance
realT clVariance(realT &varErr, realT &varNoise, const std::vector< realT > &PSDerr, const std::vector< realT > &PSDnoise, realT g)
Calculate the closed loop variance given open-loop PSDs and gain.
Definition clGainOpt.hpp:955
mx::AO::analysis::clGainOpt::clETF
complexT clETF(int fi, realT g)
Return the closed loop error transfer function (ETF) at frequency f for gain g.
Definition clGainOpt.hpp:844
mx::AO::analysis::clGainOpt::f
void f(realT *newF, size_t nF)
Set the vector of frequencies.
Definition clGainOpt.hpp:658
mx::AO::analysis::clGainOpt::clTF2
void clTF2(realT &ETF, realT &NTF, int fi, realT g)
Return the norm of the closed loop transfer functions at frequency f for gain g.
Definition clGainOpt.hpp:924
mx::AO::analysis::fourierTemporalPSD
Class to manage the calculation of temporal PSDs of the Fourier modes in atmospheric turbulence.
Definition fourierTemporalPSD.hpp:110
mx::AO::analysis::fourierTemporalPSD::_layer_i
int _layer_i
The index of the current layer.
Definition fourierTemporalPSD.hpp:134
mx::AO::analysis::fourierTemporalPSD::_w
gsl_integration_workspace * _w
Workspace for the gsl integrators, allocated to WSZ if constructed as worker (with allocate == true).
Definition fourierTemporalPSD.hpp:140
mx::AO::analysis::fourierTemporalPSD::analyzePSDGrid
int analyzePSDGrid(const std::string &subDir, const std::string &psdDir, int mnMax, int mnCon, realT gfixed, int lpNc, realT lpRegPrecision, std::vector< realT > &mags, int lifetimeTrials=0, bool ucLifeTs=false, bool writePSDs=false, bool writeXfer=false)
Analyze a PSD grid under closed-loop control.
Definition fourierTemporalPSD.hpp:759
mx::AO::analysis::fourierTemporalPSD::multiLayerPSD
int multiLayerPSD(std::vector< realT > &PSD, std::vector< realT > &freq, realT m, realT n, int p, realT fmax=0)
Calculate the temporal PSD for a Fourier mode in a multi-layer model.
Definition fourierTemporalPSD.hpp:655
mx::AO::analysis::fourierTemporalPSD::getGridFreq
int getGridFreq(std::vector< realT > &freq, const std::string &dir)
Get the frequency scale for a PSD grid.
Definition fourierTemporalPSD.hpp:2196
mx::AO::analysis::fourierTemporalPSD::m_f0
realT m_f0
the Berdja boiling parameter
Definition fourierTemporalPSD.hpp:130
mx::AO::analysis::fourierTemporalPSD::getGridPSD
int getGridPSD(std::vector< realT > &psd, const std::string &dir, int m, int n)
Get a single PSD from a PSD grid.
Definition fourierTemporalPSD.hpp:2204
mx::AO::analysis::fourierTemporalPSD::relTol
realT relTol()
Get the current relative tolerance.
Definition fourierTemporalPSD.hpp:465
mx::AO::analysis::fourierTemporalPSD::_absTol
realT _absTol
The absolute tolerance to use in the GSL integrator.
Definition fourierTemporalPSD.hpp:142
mx::AO::analysis::fourierTemporalPSD::m_modeCoeffs
Eigen::Array< realT, -1, -1 > m_modeCoeffs
Coeeficients of the projection onto the Fourier modes.
Definition fourierTemporalPSD.hpp:147
mx::AO::analysis::fourierTemporalPSD::m_cq
realT m_cq
The cosine of the wind direction.
Definition fourierTemporalPSD.hpp:123
mx::AO::analysis::fourierTemporalPSD::m_aosys
aosysT * m_aosys
Pointer to an AO system structure.
Definition fourierTemporalPSD.hpp:118
mx::AO::analysis::fourierTemporalPSD::m_spatialFilter
realT m_spatialFilter
Flag indicating if a spatial filter is applied.
Definition fourierTemporalPSD.hpp:125
mx::AO::analysis::fourierTemporalPSD::_relTol
realT _relTol
The relative tolerance to use in the GSL integrator.
Definition fourierTemporalPSD.hpp:143
mx::AO::analysis::fourierTemporalPSD::realT
_realT realT
The type for arithmetic.
Definition fourierTemporalPSD.hpp:112
mx::AO::analysis::fourierTemporalPSD::complexT
std::complex< realT > complexT
The complex type for arithmetic.
Definition fourierTemporalPSD.hpp:115
mx::AO::analysis::fourierTemporalPSD::singleLayerPSD
int singleLayerPSD(std::vector< realT > &PSD, std::vector< realT > &freq, realT m, realT n, int layer_i, int p, realT fmax=0)
Calculate the temporal PSD for a Fourier mode for a single layer.
Definition fourierTemporalPSD.hpp:494
mx::AO::analysis::fourierTemporalPSD::m_sq
realT m_sq
The sine of the wind direction.
Definition fourierTemporalPSD.hpp:124
mx::AO::analysis::fourierTemporalPSD::makePSDGrid
void makePSDGrid(const std::string &dir, int mnMax, realT dFreq, realT maxFreq, realT fmax=0)
Calculate PSDs over a grid of spatial frequencies.
Definition fourierTemporalPSD.hpp:671
mx::AO::analysis::fourierTemporalPSD::initialize
void initialize()
Initialize parameters to default values.
Definition fourierTemporalPSD.hpp:437
mx::AO::analysis::fourierTemporalPSD::fastestPeak
realT fastestPeak(int m, int n)
Determine the frequency of the highest V-dot-k peak.
Definition fourierTemporalPSD.hpp:471
mx::AO::analysis::fourierTemporalPSD::m_f
realT m_f
the current temporal frequency
Definition fourierTemporalPSD.hpp:120
mx::AO::analysis::fourierTemporalPSD::intensityPSD
int intensityPSD(const std::string &subDir, const std::string &psdDir, const std::string &CvdPath, int mnMax, int mnCon, std::vector< realT > &mags, int lifetimeTrials, bool writePSDs)
Definition fourierTemporalPSD.hpp:1510
mx::AO::analysis::fourierTemporalPSD::m_m
realT m_m
the spatial frequency m index
Definition fourierTemporalPSD.hpp:121
mx::AO::analysis::fourierTemporalPSD::m_mode_i
int m_mode_i
Projected basis mode index.
Definition fourierTemporalPSD.hpp:145
mx::AO::analysis::fourierTemporalPSD::m_n
realT m_n
the spatial frequency n index
Definition fourierTemporalPSD.hpp:122
mx::AO::analysis::fourierTemporalPSD::absTol
realT absTol()
Get the current absolute tolerance.
Definition fourierTemporalPSD.hpp:453
mx::sigproc::psdVarMean
Calculate the variance of the mean for a process given its PSD.
Definition psdVarMean.hpp:161
varmapToImage.hpp
A utility to convert a wavefront variance map to an intensity image.
wfsNoisePSD.hpp
Declares and defines a function to calculate the measurement noise PSD.