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 "../../mxlib.hpp"
41#include "../../math/constants.hpp"
42#include "../../math/func/jinc.hpp"
43#include "../../math/func/airyPattern.hpp"
44#include "../../math/vectorUtils.hpp"
45#include "../../ioutils/fits/fitsFile.hpp"
46#include "../../sigproc/fourierModes.hpp"
47#include "../../sigproc/psdVarMean.hpp"
48#include "../../ioutils/stringUtils.hpp"
49#include "../../ioutils/readColumns.hpp"
50#include "../../ioutils/binVector.hpp"
51#include "../../ioutils/fileUtils.hpp"
52
53#include "../../ipc/ompLoopWatcher.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 internal::mxlib_error_report(error_t::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 = std::format("{}/psd_{}_{}.binv",psddir,m,n);
749 // psddir + '/' + "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( 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 = std::format("{}/outputPSDS_{}_si",dir, mags[s]);
847 //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si";
848 mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
849
850 if( doLP )
851 {
852 std::string psdOutDir = std::format("{}/outputPSDS_{}_lp",dir, mags[s]);
853 //dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp";
854 mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
855 }
856 }
857 }
858
859 m_aosys->beta_p( 1, 1 );
860
861 ipc::ompLoopWatcher<> watcher( nModes * mags.size(), std::cout );
862
863 for( size_t s = 0; s < mags.size(); ++s )
864 {
865 m_aosys->starMag( mags[s] );
866
867 // In non-parallel space calculate OG=Strehl
868
869 realT opticalGain{ 1.0 };
870
871 if( m_strehlOG )
872 {
873 // Iterative optical gain
874 /// \todo need upstream NCP and CP NCP and NCP NCP
875 realT ncp = m_aosys->ncp_wfe(); // save ncp and then set it to zero for this part.
876 m_aosys->ncp_wfe( 0 );
877
878 realT lam_sci = m_aosys->lam_sci();
879 m_aosys->lam_sci( m_aosys->lam_wfs() );
880
881 realT S = m_aosys->strehl();
882 std::cerr << S << "\n";
883
884 for( int s = 0; s < 4; ++s )
885 {
886 m_aosys->opticalGain( S );
887 m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
888 S = m_aosys->strehl();
889 std::cerr << S << "\n";
890 }
891
892 opticalGain = S;
893
894 m_aosys->lam_sci( lam_sci );
895 m_aosys->ncp_wfe( ncp );
896 }
897
898 m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
899 realT strehl = m_aosys->strehl();
900
901#pragma omp parallel
902 {
903 realT localMag = mags[s];
904
905 realT var0;
906
907 realT gopt, var;
908
909 realT gopt_lp, var_lp;
910
911 std::vector<realT> tfreq; // The frequency scale of the PSDs
912 std::vector<realT> tPSDp; // The open-loop turbulence PSD for a Fourier mode
913 std::vector<realT>
914 tPSDpPOL; // The pseudo-open-loop turbulence PSD for a Fourier mode, with optical gain effects included
915
916 std::vector<realT> tfreqHF; // The above-Nyquist frequencies, saved if outputing the PSDS.
917 std::vector<realT> tPSDpHF; // The above-Nyquist component of the open-loop PSD, saved if outputing the
918 // PSDs.
919
920 //**< Get the frequency grid, and nyquist limit it to f_s/2
921 getGridPSD( tfreq, tPSDp, psdDir, 0, 1 ); // To get the freq grid
922
923 size_t imax = 0;
924 while( tfreq[imax] <= 0.5 * fs )
925 {
926 ++imax;
927 if( imax > tfreq.size() - 1 )
928 break;
929 }
930
931 if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
932 {
933 ++imax;
934 }
935
936 if( writePSDs )
937 {
938 tfreqHF.assign( tfreq.begin(), tfreq.end() );
939 }
940
941 tfreq.erase( tfreq.begin() + imax, tfreq.end() );
942
943 tPSDpPOL.resize( tfreq.size() ); // pre=allocate
944 //**>
945
946 std::vector<realT> tPSDn; // The open-loop WFS noise PSD
947 tPSDn.resize( tfreq.size() );
948
949 //**< Setup the controllers
950 mx::AO::analysis::clAOLinearPredictor<realT> tflp;
951 tflp.m_precision0 = lpRegPrecision;
952 /*tflp.m_min_sc0 = 0;
953 tflp.m_max_sc0 = 1000;*/
954
955 mx::AO::analysis::clGainOpt<realT> go_si( tauWFS, deltaTau );
956 mx::AO::analysis::clGainOpt<realT> go_lp( tauWFS, deltaTau );
957
958 go_si.f( tfreq );
959 go_lp.f( tfreq );
960
961 realT gmax = 0;
962 realT gmax_lp = 0;
963 //**>
964
965 int m, n;
966
967 //**< For use in lifetime calculations
968 sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
969 std::vector<std::complex<realT>> ETFxn;
970 std::vector<std::complex<realT>> NTFxn;
971
972 if( lifetimeTrials > 0 )
973 {
974 ETFxn.resize( tfreq.size() );
975 NTFxn.resize( tfreq.size() );
976
977 if( writeXfer )
978 {
979 std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]);
980 //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
981 ioutils::createDirectories( tfOutFile );
982 }
983
984 if( doLP )
985 {
986 if( writeXfer )
987 {
988 std::string tfOutFile = std::format("{}/outputTF_{}_lp", dir, mags[s]);
989 //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
990 ioutils::createDirectories( tfOutFile );
991 }
992 }
993 }
994
995//**>
996
997/*#pragma omp critical
998std::cerr << __FILE__ << " " << __LINE__ << "\n";
999*/
1000// want to schedule dynamic with small chunks so maximal processor usage,
1001// otherwise we can end up with a small number of cores being used at the end
1002#pragma omp for schedule( dynamic, 5 )
1003 for( size_t i = 0; i < nModes; ++i )
1004 {
1005 // Determine the spatial frequency at this step
1006 m = fms[2 * i].m;
1007 n = fms[2 * i].n;
1008
1009 if( fabs( (realT)m / m_aosys->D() ) >= m_aosys->spatialFilter_ku() ||
1010 fabs( (realT)n / m_aosys->D() ) >= m_aosys->spatialFilter_kv() )
1011 {
1012 gains( mnMax + m, mnMax + n ) = 0;
1013 gains( mnMax - m, mnMax - n ) = 0;
1014
1015 gains_lp( mnMax + m, mnMax + n ) = 0;
1016 gains_lp( mnMax - m, mnMax - n ) = 0;
1017
1018 vars( mnMax + m, mnMax + n ) = 0;
1019 vars( mnMax - m, mnMax - n ) = 0;
1020
1021 vars_lp( mnMax + m, mnMax + n ) = 0;
1022 vars_lp( mnMax - m, mnMax - n ) = 0;
1023 speckleLifetimes( mnMax + m, mnMax + n ) = 0;
1024 speckleLifetimes( mnMax - m, mnMax - n ) = 0;
1025 speckleLifetimes_lp( mnMax + m, mnMax + n ) = 0;
1026 speckleLifetimes_lp( mnMax - m, mnMax - n ) = 0;
1027 }
1028 else
1029 {
1030
1031 realT k = sqrt( m * m + n * n ) / m_aosys->D();
1032
1033 //**< Get the open-loop turb. PSD
1034 getGridPSD( tPSDp, psdDir, m, n );
1035
1036 // Get integral of entire open-loop PSD
1037 var0 = sigproc::psdVar( tfreq, tPSDp );
1038
1039 if( writePSDs )
1040 {
1041 tPSDpHF.assign( tPSDp.begin() + imax, tPSDp.end() );
1042 }
1043
1044 // erase points above Nyquist limit
1045 tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() );
1046
1047 // And now determine the variance which has been erased.
1048 // limVar is the out-of-band variance, which we add back in for completeness
1049 realT limVar = 0; // var0 - sigproc::psdVar( tfreq, tPSDp);
1050
1051 // And construct the POL psd
1052 if( m_uncorrectedOG )
1053 {
1054 for( size_t n = 0; n < tPSDp.size(); ++n )
1055 {
1056 tPSDpPOL[n] = tPSDp[n] * pow( opticalGain, 2 );
1057 }
1058 }
1059 else
1060 {
1061 for( size_t n = 0; n < tPSDp.size(); ++n )
1062 {
1063 tPSDpPOL[n] = tPSDp[n];
1064 }
1065 }
1066 //**>
1067
1068 //**< Determine if we're inside the hardwarecontrol limit
1069 bool inside = false;
1070
1071 if( m_aosys->circularLimit() )
1072 {
1073 if( m * m + n * n <= mnCon * mnCon )
1074 inside = true;
1075 }
1076 else
1077 {
1078 if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
1079 inside = true;
1080 }
1081 //**>
1082
1083 /* This is to select specific points for troubleshooting*/
1084 // if( !( ( m ==-2 && n == 84 ) || (m==-2 && n == 2800)) ) inside = false;
1085
1086 if( inside )
1087 {
1088 // Get the WFS noise PSD (which is already resized to match tfreq)
1089 wfsNoisePSD<realT>( tPSDn,
1090 m_aosys->beta_p( m, n ) / sqrt( opticalGain ),
1091 m_aosys->Fg( localMag ),
1092 tauWFS,
1093 m_aosys->npix_wfs( (size_t)0 ),
1094 m_aosys->Fbg( (size_t)0 ),
1095 m_aosys->ron_wfs( (size_t)0 ) );
1096
1097 gmax = 0;
1098 if( gfixed > 0 )
1099 {
1100 gopt = gfixed;
1101 var = go_si.clVariance( tPSDp, tPSDn, gopt );
1102 }
1103 else
1104 {
1105 // Calculate gain using the POL PSD
1106 gopt = go_si.optGainOpenLoop( var, tPSDpPOL, tPSDn, gmax, true );
1107
1108 if( m_uncorrectedOG )
1109 {
1110 gopt *= opticalGain;
1111 }
1112
1113 // But use the variance from the actual POL
1114 var = go_si.clVariance( tPSDp, tPSDn, gopt );
1115
1116 // Output gain curve for this mode (for troubleshooting)
1117 /*#pragma omp critical
1118 {
1119 std::string foutn = "gcurve_";
1120 foutn += std::to_string(m) + "_" + std::to_string(n) + ".dat";
1121
1122 std::ofstream foutf(foutn);
1123
1124 for(size_t n = 0; n < 10000; ++n)
1125 {
1126 realT gg = (1.0*n)/10000.;
1127 foutf << gg << " " << go_si.clVariance(tPSDp, tPSDn, gg) << "\n";
1128 }
1129
1130 foutf.close();
1131
1132 std::cerr << "\n" << gmax << " " << gopt << " " << var << " " << go_si.clVariance(tPSDp,
1133 tPSDn, 0.64) << "\n";
1134 }*/
1135 }
1136
1137 var += limVar;
1138
1139 if( doLP )
1140 {
1141 realT min_sc;
1142 int rv = tflp.regularizeCoefficients( gmax_lp,
1143 gopt_lp,
1144 var_lp,
1145 min_sc,
1146 go_lp,
1147 tPSDpPOL,
1148 tPSDn,
1149 lpNc );
1150
1151 if( rv < 0 )
1152 {
1153 std::cerr
1154 << "fourierTemporalPSD::analyzePSDGrid: regularizeCoefficients returned error ";
1155 std::cerr << rv << ' ';
1156 std::cerr << __FILE__ << ' ' << __LINE__ << '\n';
1157 }
1158
1159 for( int n = 0; n < lpNc; ++n )
1160 {
1161 lpC( i, n ) = go_lp.a( n );
1162 }
1163
1164 if( m_uncorrectedOG )
1165 {
1166 go_lp.aScale( opticalGain );
1167 go_lp.bScale( opticalGain );
1168 gopt_lp *= opticalGain;
1169 }
1170
1171 var_lp = go_lp.clVariance( tPSDp, tPSDn, gopt_lp );
1172 var_lp += limVar;
1173 }
1174 else
1175 {
1176 gopt_lp = 0;
1177 }
1178 }
1179 else
1180 {
1181 // Zero the noise PSD
1182 tPSDn.assign( tPSDn.size(), 0.0 );
1183
1184 gopt = 0;
1185 var = var0;
1186 var_lp = var0;
1187 gopt_lp = 0;
1188 go_lp.a( std::vector<realT>( { 1 } ) );
1189 go_lp.b( std::vector<realT>( { 1 } ) );
1190 }
1191
1192 //**< Determine if closed-loop is making a difference:
1193
1194 if( gopt > 0 && var > var0 )
1195 {
1196 gopt = 0;
1197 var = var0;
1198 }
1199
1200 if( gopt_lp > gopt && var_lp > var )
1201 {
1202 //Set LP to SI (or off if SI is off)
1203 gopt_lp = gopt;
1204 var_lp = var;
1205 go_lp.a( std::vector<realT>( { 1 } ) );
1206 go_lp.b( std::vector<realT>( { 1 } ) );
1207 }
1208 //**>
1209
1210 //**< Fill in the gain and variance maps
1211 gains( mnMax + m, mnMax + n ) = gopt;
1212 gains( mnMax - m, mnMax - n ) = gopt;
1213
1214 gains_lp( mnMax + m, mnMax + n ) = gopt_lp;
1215 gains_lp( mnMax - m, mnMax - n ) = gopt_lp;
1216
1217 vars( mnMax + m, mnMax + n ) = var;
1218 vars( mnMax - m, mnMax - n ) = var;
1219
1220 vars_lp( mnMax + m, mnMax + n ) = var_lp;
1221 vars_lp( mnMax - m, mnMax - n ) = var_lp;
1222 //**>
1223
1224 //**< Calculate Speckle Lifetimes
1225 if( ( lifetimeTrials > 0 || writeXfer ) && ( uncontrolledLifetimes || inside ) )
1226 {
1227 std::vector<realT> spfreq, sppsd;
1228
1229 if( gopt > 0 )
1230 {
1231 for( size_t i = 0; i < tfreq.size(); ++i )
1232 {
1233 ETFxn[i] = go_si.clETF( i, gopt );
1234 NTFxn[i] = go_si.clNTF( i, gopt );
1235 }
1236 }
1237 else
1238 {
1239 for( size_t i = 0; i < tfreq.size(); ++i )
1240 {
1241 ETFxn[i] = 1;
1242 NTFxn[i] = 0;
1243 }
1244 }
1245
1246 if( writeXfer )
1247 {
1248 std::string tfOutFile = std::format("{}/outputTF_{}_si", dir, mags[s]);
1249 // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
1250
1251 std::string etfOutFile = std::format("{}/etf_{}_{}.binv", tfOutFile,m, n);
1252 //tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' +
1253 // ioutils::convert ToString( n ) + ".binv";
1254 ioutils::writeBinVector( etfOutFile, ETFxn );
1255
1256 std::string ntfOutFile = std::format("{}/ntf_{}_{}.binv",tfOutFile, m, n);
1257 //tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' +
1258 // ioutils::convert ToString( n ) + ".binv";
1259 ioutils::writeBinVector( ntfOutFile, NTFxn );
1260
1261 if( i == 0 ) // Write freq on the first one
1262 {
1263 std::string fOutFile = tfOutFile + "freq.binv";
1264 ioutils::writeBinVector( fOutFile, tfreq );
1265 }
1266 }
1267
1268 if( lifetimeTrials > 0 )
1269 {
1270 speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1271 realT spvar = sigproc::psdVar( spfreq, sppsd );
1272
1273 realT splifeT = 100.0;
1274 realT error;
1275
1276 realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1277
1278 speckleLifetimes( mnMax + m, mnMax + n ) = tau;
1279 speckleLifetimes( mnMax - m, mnMax - n ) = tau;
1280 }
1281
1282 if( doLP )
1283 {
1284 if( gopt_lp > 0 )
1285 {
1286 for( size_t i = 0; i < tfreq.size(); ++i )
1287 {
1288 ETFxn[i] = go_lp.clETF( i, gopt_lp );
1289 NTFxn[i] = go_lp.clNTF( i, gopt_lp );
1290 }
1291 }
1292 else
1293 {
1294 for( size_t i = 0; i < tfreq.size(); ++i )
1295 {
1296 ETFxn[i] = 1;
1297 NTFxn[i] = 0;
1298 }
1299 }
1300
1301 if( writeXfer )
1302 {
1303 std::string tfOutFile = std::format("{}/outputTF_{}_lp",dir,mags[s]);
1304 //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
1305
1306 std::string etfOutFile = std::format("{}/etf_{}_{}.binv", tfOutFile, m, n);
1307 //tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' +
1308 // ioutils::convert ToString( n ) + ".binv";
1309 ioutils::writeBinVector( etfOutFile, ETFxn );
1310
1311 std::string ntfOutFile = std::format("{}/ntf_{}_{}.binv", tfOutFile, m, n);
1312 //tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' +
1313 // ioutils::convert ToString( n ) + ".binv";
1314 ioutils::writeBinVector( ntfOutFile, NTFxn );
1315
1316 if( i == 0 ) // Write freq on the first one
1317 {
1318 std::string fOutFile = tfOutFile + "freq.binv";
1319 ioutils::writeBinVector( fOutFile, tfreq );
1320 }
1321 }
1322
1323 if( lifetimeTrials > 0 )
1324 {
1325 speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1326 realT spvar = sigproc::psdVar( spfreq, sppsd );
1327
1328 realT splifeT = 100.0;
1329 realT error;
1330
1331 realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1332
1333 speckleLifetimes_lp( mnMax + m, mnMax + n ) = tau;
1334 speckleLifetimes_lp( mnMax - m, mnMax - n ) = tau;
1335 }
1336 } // if(doLP)
1337 } // if( (lifetimeTrials > 0 || writeXfer) && ( uncontrolledLifetimes || inside ))
1338 //**>
1339
1340 // Calculate the controlled PSDs and output
1341 if( writePSDs )
1342 {
1343 std::string psdOutFile = std::format("{}/outputPSDs_{}_si/psd_{}_{}.binv",dir,mags[s],m,n);
1344
1345 // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/";
1346 //psdOutFile +=
1347 // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1348
1349 std::vector<realT> psdOut( tPSDp.size() + tPSDpHF.size() );
1350
1351 // Calculate the output PSD if gains are applied
1352 if( gopt > 0 )
1353 {
1354 realT ETF, NTF;
1355
1356 for( size_t i = 0; i < tfreq.size(); ++i )
1357 {
1358 go_si.clTF2( ETF, NTF, i, gopt );
1359 psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
1360 }
1361
1362 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1363 {
1364 psdOut[tfreq.size() + i] = tPSDpHF[i];
1365 }
1366 }
1367 else // otherwise just copy
1368 {
1369 for( size_t i = 0; i < tfreq.size(); ++i )
1370 {
1371 psdOut[i] = tPSDp[i];
1372 }
1373
1374 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1375 {
1376 psdOut[tfreq.size() + i] = tPSDpHF[i];
1377 }
1378 }
1379
1380 ioutils::writeBinVector( psdOutFile, psdOut );
1381
1382 if( i == 0 ) // Write freq on the first one
1383 {
1384 psdOutFile = std::format("{}/outputPSDs_{}_si/freq.binv", dir, mags[s]);
1385 // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/freq.binv";
1386 //ioutils::writeBinVector( psdOutFile, tfreqHF );
1387 }
1388
1389 if( doLP )
1390 {
1391 std::string psdOutFile = std::format("{}/outputPSDs_{}_lp/psd_{}_{}.binv",dir,mags[s],m,n);
1392
1393 //psdOutFile = dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/";
1394 //psdOutFile +=
1395 // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1396
1397 // Calculate the output PSD if gains are applied
1398 if( gopt_lp > 0 )
1399 {
1400 realT ETF, NTF;
1401
1402 for( size_t i = 0; i < tfreq.size(); ++i )
1403 {
1404 go_lp.clTF2( ETF, NTF, i, gopt_lp );
1405 psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
1406 }
1407 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1408 {
1409 psdOut[tfreq.size() + i] = tPSDpHF[i];
1410 }
1411 }
1412 else // otherwise just copy
1413 {
1414 for( size_t i = 0; i < tfreq.size(); ++i )
1415 {
1416 psdOut[i] = tPSDp[i];
1417 }
1418
1419 for( size_t i = 0; i < tPSDpHF.size(); ++i )
1420 {
1421 psdOut[tfreq.size() + i] = tPSDpHF[i];
1422 }
1423 }
1424
1425 ioutils::writeBinVector( psdOutFile, psdOut );
1426
1427 if( i == 0 )
1428 {
1429 psdOutFile = std::format("{}/outputPSDs_{}_lp/freq.binv", dir, mags[s]);
1430 //psdOutFile =
1431 // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/freq.binv";
1432 //ioutils::writeBinVector( psdOutFile, tfreq );
1433 }
1434 }
1435 }
1436 }
1437 watcher.incrementAndOutputStatus();
1438
1439 } // omp for i..nModes
1440 } // omp Parallel
1441
1442 Eigen::Array<realT, -1, -1> cim;
1443
1444 fits::fitsFile<realT> ff;
1445 std::string fn = std::format("{}/gainmap_{}_si.fits",dir, mags[s]);
1446 //dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
1447 ff.write( fn, gains );
1448
1449 fn = std::format("{}/varmap_{}_si.fits",dir, mags[s]);
1450 //dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
1451 ff.write( fn, vars );
1452
1453 cim = vars;
1454
1455 realT Ssi = exp( -1 * cim.sum() );
1456 S_si.push_back( strehl );
1457 cim /= strehl;
1458
1459 fn = std::format("{}/contrast_{}_si.fits",dir, mags[s]);
1460 //dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
1461 ff.write( fn, cim );
1462
1463 if( lifetimeTrials > 0 )
1464 {
1465 fn = std::format("{}/speckleLifetimes_{}_si.fits",dir, mags[s]);
1466 //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
1467 ff.write( fn, speckleLifetimes );
1468 }
1469
1470 if( doLP )
1471 {
1472 fn = std::format("{}/gainmap_{}_lp.fits",dir, mags[s]);
1473 //dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
1474 ff.write( fn, gains_lp );
1475
1476 fn = std::format("{}/lpcmap_{}_lp.fits",dir, mags[s]);
1477 //dir + "/lpcmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
1478 ff.write( fn, lpC );
1479
1480 fn = std::format("{}/varmap_{}_lp.fits",dir, mags[s]);
1481 //dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
1482 ff.write( fn, vars_lp );
1483
1484 cim = vars_lp;
1485
1486 // Scale Strehl by the ratio of total variance
1487 realT Slp = strehl * exp( -1 * cim.sum() ) /
1488 Ssi; // This is a hack until we do a real fitting error calculation or something
1489 S_lp.push_back( Slp );
1490 cim /= Slp;
1491
1492 fn = std::format("{}/contrast_{}_lp.fits",dir, mags[s]);
1493 //dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
1494 ff.write( fn, cim );
1495
1496 if( lifetimeTrials > 0 )
1497 {
1498 fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]);
1499 //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
1500 ff.write( fn, speckleLifetimes_lp );
1501 }
1502 }
1503
1504 } // s (mag)
1505
1506 fn = dir + "/strehl_si.txt";
1507 fout.open( fn );
1508 for( size_t i = 0; i < mags.size(); ++i )
1509 {
1510 fout << mags[i] << " " << S_si[i] << "\n";
1511 }
1512
1513 fout.close();
1514
1515 if( doLP )
1516 {
1517 fn = dir + "/strehl_lp.txt";
1518 fout.open( fn );
1519 for( size_t i = 0; i < mags.size(); ++i )
1520 {
1521 fout << mags[i] << " " << S_lp[i] << "\n";
1522 }
1523
1524 fout.close();
1525 }
1526
1527 return 0;
1528}
1529
1530template <typename realT, typename aosysT>
1531int fourierTemporalPSD<realT, aosysT>::intensityPSD(
1532 const std::string &subDir, // sub-directory of psdDir which contains the controlled system results,
1533 // and where the lifetimes will be written.
1534 const std::string &psdDir, // directory containing the PSDS
1535 const std::string &CvdPath, // path to the covariance decomposition
1536 int mnMax,
1537 int mnCon,
1538 std::vector<realT> &mags,
1539 int lifetimeTrials,
1540 bool writePSDs )
1541{
1542
1543 std::string dir = psdDir + "/" + subDir;
1544
1545 /*** Dump Params to file ***/
1546 mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
1547
1548 std::ofstream fout;
1549 std::string fn = dir + '/' + "splife_params.txt";
1550 fout.open( fn );
1551
1552 fout << "#---------------------------\n";
1553 m_aosys->dumpAOSystem( fout );
1554 fout << "#---------------------------\n";
1555 fout << "# Analysis Parameters\n";
1556 fout << "# mnMax = " << mnMax << '\n';
1557 fout << "# mnCon = " << mnCon << '\n';
1558 fout << "# mags = ";
1559 for( size_t i = 0; i < mags.size() - 1; ++i )
1560 fout << mags[i] << ",";
1561 fout << mags[mags.size() - 1] << '\n';
1562 fout << "# lifetimeTrials = " << lifetimeTrials << '\n';
1563 // fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n';
1564 fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n';
1565
1566 fout.close();
1567
1568 realT fs = 1.0 / m_aosys->tauWFS();
1569 realT tauWFS = m_aosys->tauWFS();
1570 realT deltaTau = m_aosys->deltaTau();
1571
1572 //** Get the Fourier Grid
1573 std::vector<sigproc::fourierModeDef> fms;
1574
1575 sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax );
1576 size_t nModes = 0.5 * fms.size();
1577 std::cerr << "nModes: " << nModes << " (" << fms.size() << ")\n";
1578
1579 Eigen::Array<realT, -1, -1> speckleLifetimes;
1580 Eigen::Array<realT, -1, -1> speckleLifetimes_lp;
1581
1582 speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1583 speckleLifetimes( mnMax, mnMax ) = 0;
1584
1585 speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1586 speckleLifetimes_lp( mnMax, mnMax ) = 0;
1587
1588 /*********************************************************************/
1589 // 0) Get the frequency grid, and nyquist limit it to f_s/2
1590 /*********************************************************************/
1591
1592 std::vector<realT> tfreq;
1593 std::vector<realT> tPSDp; // The open-loop OPD PSD
1594 std::vector<realT> tPSDn; // The open-loop WFS noise PSD
1595 std::vector<complexT> tETF;
1596 std::vector<complexT> tNTF;
1597
1598 if( getGridFreq( tfreq, psdDir ) < 0 )
1599 return -1;
1600
1601 size_t imax = 0;
1602 while( tfreq[imax] <= 0.5 * fs )
1603 {
1604 ++imax;
1605 if( imax > tfreq.size() - 1 )
1606 break;
1607 }
1608
1609 if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
1610 ++imax;
1611
1612 tfreq.erase( tfreq.begin() + imax, tfreq.end() );
1613
1614 // Now allocate memory
1615 tPSDn.resize( tfreq.size() );
1616 std::vector<std::vector<realT>> sqrtOPDPSD;
1617 sqrtOPDPSD.resize( nModes );
1618
1619 std::vector<std::vector<realT>> opdPSD;
1620 opdPSD.resize( nModes );
1621
1622 std::vector<realT> psd2sided;
1623
1624 // Store the mode variance for later normalization
1625 std::vector<realT> modeVar;
1626 modeVar.resize( nModes );
1627
1628 /*********************************************************************/
1629 // 1) Read in each PSD, and load it into the array in FFT order
1630 /*********************************************************************/
1631
1632 for( size_t i = 0; i < nModes; ++i )
1633 {
1634 // Determine the spatial frequency at this step
1635 int m = fms[2 * i].m;
1636 int n = fms[2 * i].n;
1637
1638 //**< Get the open-loop turb. PSD
1639 if( getGridPSD( tPSDp, psdDir, m, n ) < 0 )
1640 return -1;
1641 tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() ); // Nyquist limit
1642 modeVar[i] = sigproc::psdVar( tfreq, tPSDp );
1643
1644 // And now normalize
1645 sigproc::normPSD( tPSDp, tfreq, 1.0 ); // Normalize
1646 sigproc::augment1SidedPSD( psd2sided, tPSDp, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
1647
1648 opdPSD[i].resize( psd2sided.size() );
1649 sqrtOPDPSD[i].resize( psd2sided.size() );
1650
1651 for( size_t j = 0; j < psd2sided.size(); ++j )
1652 {
1653 opdPSD[i][j] = psd2sided[j] * modeVar[i];
1654 sqrtOPDPSD[i][j] = sqrt( psd2sided[j] ); // Store the square-root for later
1655 }
1656 }
1657
1658 size_t sz2Sided = psd2sided.size();
1659
1660 std::vector<realT> freq2sided;
1661 freq2sided.resize( sz2Sided );
1662 sigproc::augment1SidedPSDFreq( freq2sided, tfreq );
1663
1664 tPSDp.resize( tfreq.size() );
1665 tETF.resize( tfreq.size() );
1666 tNTF.resize( tfreq.size() );
1667
1668 std::vector<std::vector<realT>> sqrtNPSD;
1669 sqrtNPSD.resize( nModes );
1670
1671 std::vector<realT> noiseVar;
1672 noiseVar.resize( nModes );
1673
1674 std::vector<std::vector<complexT>> ETFsi;
1675 std::vector<std::vector<complexT>> ETFlp;
1676 ETFsi.resize( nModes );
1677 ETFlp.resize( nModes );
1678
1679 std::vector<std::vector<complexT>> NTFsi;
1680 std::vector<std::vector<complexT>> NTFlp;
1681 NTFsi.resize( nModes );
1682 NTFlp.resize( nModes );
1683
1684 std::string tfInFile;
1685 std::string etfInFile;
1686 std::string ntfInFile;
1687
1688 improc::eigenImage<realT> Cvd;
1689 fits::fitsFile<realT> ff;
1690 ff.read( Cvd, CvdPath );
1691
1692 std::vector<std::complex<realT>> tPSDc, psd2sidedc;
1693
1694 /*********************************************************************/
1695 // 2) Analyze each star magnitude
1696 /*********************************************************************/
1697 ipc::ompLoopWatcher<> watcher( lifetimeTrials * mags.size(), std::cout );
1698 for( size_t s = 0; s < mags.size(); ++s )
1699 {
1700 /*********************************************************************/
1701 // 2.0) Read in the transfer functions for each mode
1702 /*********************************************************************/
1703 for( size_t i = 0; i < nModes; ++i )
1704 {
1705 // Determine the spatial frequency at this step
1706 int m = fms[2 * i].m;
1707 int n = fms[2 * i].n;
1708
1709 //**< Determine if we're inside the hardwarecontrol limit
1710 bool inside = false;
1711
1712 if( m_aosys->circularLimit() )
1713 {
1714 if( m * m + n * n <= mnCon * mnCon )
1715 inside = true;
1716 }
1717 else
1718 {
1719 if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
1720 inside = true;
1721 }
1722
1723 // Get the WFS noise PSD (which is already resized to match tfreq)
1724 wfsNoisePSD<realT>( tPSDn,
1725 m_aosys->beta_p( m, n ),
1726 m_aosys->Fg( mags[s] ),
1727 tauWFS,
1728 m_aosys->npix_wfs( (size_t)0 ),
1729 m_aosys->Fbg( (size_t)0 ),
1730 m_aosys->ron_wfs( (size_t)0 ) );
1731 sigproc::augment1SidedPSD( psd2sided, tPSDn, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
1732
1733 // Pre-calculate the variance of the noise for later use
1734 noiseVar[i] = sigproc::psdVar( tfreq, tPSDn );
1735
1736 sqrtNPSD[i].resize( psd2sided.size() );
1737 for( size_t j = 0; j < psd2sided.size(); ++j )
1738 sqrtNPSD[i][j] = sqrt( psd2sided[j] );
1739
1740 ETFsi[i].resize( sz2Sided );
1741 ETFlp[i].resize( sz2Sided );
1742 NTFsi[i].resize( sz2Sided );
1743 NTFlp[i].resize( sz2Sided );
1744
1745 if( inside )
1746 {
1747 tfInFile = std::format("{}/outputTF_{}_si",dir, mags[s]);
1748 //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
1749
1750 etfInFile = std::format("{}/etf_{}_{}.binv",tfInFile, m, n);
1751 //tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1752 ioutils::readBinVector( tPSDc, etfInFile );
1753 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1754 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1755 ETFsi[i][j] = psd2sidedc[j];
1756
1757 ntfInFile = std::format("{}/ntf_{}_{}.binv",tfInFile, m, n);
1758 //tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1759 ioutils::readBinVector( tPSDc, ntfInFile );
1760 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1761 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1762 NTFsi[i][j] = psd2sidedc[j];
1763
1764 tfInFile = std::format("{}/outputTF_{}_lp",dir, mags[s]);
1765 //dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
1766
1767 etfInFile = std::format("{}/etf_{}_{}.binv",tfInFile, m, n);
1768 //tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1769 ioutils::readBinVector( tPSDc, etfInFile );
1770 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1771 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1772 ETFlp[i][j] = psd2sidedc[j];
1773
1774 ntfInFile = std::format("{}/ntf_{}_{}.binv",tfInFile, m, n);
1775 //tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1776 ioutils::readBinVector( tPSDc, ntfInFile );
1777 sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
1778 for( size_t j = 0; j < psd2sidedc.size(); ++j )
1779 NTFlp[i][j] = psd2sidedc[j];
1780 }
1781 else
1782 {
1783 for( int q = 0; q < ETFsi.size(); ++q )
1784 {
1785 ETFsi[i][q] = 1;
1786 NTFsi[i][q] = 0;
1787 ETFlp[i][q] = 1;
1788 NTFlp[i][q] = 0;
1789 }
1790 }
1791 }
1792
1793 sigproc::averagePeriodogram<realT> tavgPgram(
1794 sz2Sided / 1.,
1795 1 / fs ); // this is just to get the size right, per-thread instances below
1796 std::vector<std::vector<realT>> spPSDs;
1797 spPSDs.resize( nModes );
1798 for( size_t pp = 0; pp < spPSDs.size(); ++pp )
1799 {
1800 spPSDs[pp].resize( tavgPgram.size() );
1801 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
1802 spPSDs[pp][nn] = 0;
1803 }
1804
1805 std::vector<std::vector<realT>> spPSDslp;
1806 spPSDslp.resize( nModes );
1807 for( size_t pp = 0; pp < spPSDslp.size(); ++pp )
1808 {
1809 spPSDslp[pp].resize( tavgPgram.size() );
1810 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
1811 spPSDslp[pp][nn] = 0;
1812 }
1813
1814#pragma omp parallel
1815 {
1816 // Normally distributed random numbers
1817 math::normDistT<realT> normVar;
1818 normVar.seed();
1819
1820 // FFTs for going to Fourier domain and back to time domain.
1821 math::ft::fftT<realT, std::complex<realT>, 1, 0> fftF( sqrtOPDPSD[0].size() );
1822 math::ft::fftT<std::complex<realT>, realT, 1, 0> fftB( sqrtOPDPSD[0].size(), math::ft::dir::backward );
1823
1824 // Working memory
1825 std::vector<std::complex<realT>> tform1( sqrtOPDPSD[0].size() );
1826 std::vector<std::complex<realT>> tform2( sqrtOPDPSD[0].size() );
1827 std::vector<std::complex<realT>> Ntform1( sqrtOPDPSD[0].size() );
1828 std::vector<std::complex<realT>> Ntform2( sqrtOPDPSD[0].size() );
1829
1830 std::vector<std::complex<realT>> tform1lp( sqrtOPDPSD[0].size() );
1831 std::vector<std::complex<realT>> tform2lp( sqrtOPDPSD[0].size() );
1832 std::vector<std::complex<realT>> Ntform1lp( sqrtOPDPSD[0].size() );
1833 std::vector<std::complex<realT>> Ntform2lp( sqrtOPDPSD[0].size() );
1834
1835 // OPD-PSD filter
1836 sigproc::psdFilter<realT, 1> pfilt;
1837 pfilt.psdSqrt( &sqrtOPDPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
1838
1839 // Noise-PSD filter
1840 sigproc::psdFilter<realT, 1> nfilt;
1841 nfilt.psdSqrt( &sqrtNPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
1842
1843 // The h time-series
1844 std::vector<std::vector<realT>> hts;
1845 hts.resize( 2 * nModes );
1846
1847 // The correlated h time-series
1848 std::vector<std::vector<realT>> htsCorr;
1849 htsCorr.resize( 2 * nModes );
1850
1851 for( size_t pp = 0; pp < hts.size(); ++pp )
1852 {
1853 hts[pp].resize( sqrtOPDPSD[0].size() );
1854 htsCorr[pp].resize( sqrtOPDPSD[0].size() );
1855 }
1856
1857 // The noise time-serieses
1858 std::vector<realT> N_n;
1859 N_n.resize( sz2Sided );
1860
1861 std::vector<realT> N_nm;
1862 N_nm.resize( sz2Sided );
1863
1864 // Periodogram averager
1865 sigproc::averagePeriodogram<realT> avgPgram( sz2Sided / 1., 1 / fs ); //, 0, 1);
1866 avgPgram.win( sigproc::window::hann );
1867
1868 // The periodogram output
1869 std::vector<realT> tpgram( avgPgram.size() );
1870
1871 // Holds the speckle time-series
1872 improc::eigenCube<realT> spTS;
1873 spTS.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
1874
1875 improc::eigenCube<realT> spTSlp;
1876 spTSlp.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
1877
1878// Here's where the big loop of n-trials should start
1879#pragma omp for
1880 for( int zz = 0; zz < lifetimeTrials; ++zz )
1881 {
1882 std::complex<realT> scale = exp( std::complex<realT>( 0, math::half_pi<realT>() ) ) /
1883 std::complex<realT>( ( tform1.size() ), 0 );
1884
1885 /*********************************************************************/
1886 // 2.1) Generate filtered noise for each mode, with temporal phase shifts at each spatial frequency
1887 /*********************************************************************/
1888 for( size_t pp = 0; pp < nModes; ++pp )
1889 {
1890 // Fill in standard normal noise
1891 for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
1892 {
1893 hts[2 * pp][nn] = normVar;
1894 }
1895
1896 // Set sqrt(PSD), just a pointer switch
1897 pfilt.psdSqrt( &sqrtOPDPSD[pp], tfreq[1] - tfreq[0] );
1898
1899 // And now filter the noise to a time-series of h
1900 pfilt( hts[2 * pp] );
1901
1902 /**/
1903 // Then construct 2nd mode with temporal shift
1904 fftF( tform1.data(), hts[2 * pp].data() );
1905
1906 // Apply the phase shift to form the 2nd time series
1907 for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
1908 tform1[nn] = tform1[nn] * scale;
1909
1910 fftB( hts[2 * pp + 1].data(), tform1.data() );
1911 /**/
1912 }
1913 //** At this point we have correlated time-series, with the correct temporal PSD, but not yet spatially
1914 // correlated
1915
1916 /*********************************************************************/
1917 // 2.2) Correlate the time-series for each mode
1918 /*********************************************************************/
1919 // #pragma omp parallel for
1920 for( size_t pp = 0; pp < hts.size(); ++pp )
1921 {
1922 for( size_t nn = 0; nn < hts[0].size(); ++nn )
1923 {
1924 htsCorr[pp][nn] = 0;
1925 }
1926
1927 for( size_t qq = 0; qq <= pp; ++qq )
1928 {
1929 realT cvd = Cvd( qq, pp );
1930 realT *d1 = htsCorr[pp].data();
1931 realT *d2 = hts[qq].data();
1932 for( size_t nn = 0; nn < hts[0].size(); ++nn )
1933 {
1934 d1[nn] += d2[nn] * cvd;
1935 }
1936 }
1937 }
1938
1939 /*
1940 for(size_t pp=0; pp < hts.size(); ++pp)
1941 {
1942 for(size_t nn=0; nn< hts[0].size(); ++nn)
1943 {
1944 htsCorr[pp][nn] = hts[pp][nn];
1945 }
1946 }*/
1947
1948 /*********************************************************************/
1949 // 2.2.a) Re-normalize b/c the correlation step above does not result in correct variances
1950 ///\todo should be able to scale the covar by r0, and possibly D
1951 /*********************************************************************/
1952 for( size_t pp = 0; pp < nModes; ++pp )
1953 {
1954 math::vectorMeanSub( htsCorr[2 * pp] );
1955 math::vectorMeanSub( htsCorr[2 * pp + 1] );
1956
1957 realT var = math::vectorVariance( htsCorr[2 * pp] );
1958 realT norm = sqrt( modeVar[pp] / var );
1959 for( size_t nn = 0; nn < htsCorr[2 * pp].size(); ++nn )
1960 htsCorr[2 * pp][nn] *= norm;
1961
1962 var = math::vectorVariance( htsCorr[2 * pp + 1] );
1963 norm = sqrt( modeVar[pp] / var );
1964 for( size_t nn = 0; nn < htsCorr[2 * pp + 1].size(); ++nn )
1965 htsCorr[2 * pp + 1][nn] *= norm;
1966 }
1967
1968 scale = std::complex<realT>( tform1.size(), 0 );
1969
1970 /*********************************************************************/
1971 // 2.3) Generate speckle intensity time-series
1972 /*********************************************************************/
1973 for( size_t pp = 0; pp < nModes; ++pp )
1974 {
1975 // Now we take them back to the FD and apply the xfer
1976 // and add the noise
1977
1978 fftF( tform1.data(), htsCorr[2 * pp].data() );
1979 fftF( tform2.data(), htsCorr[2 * pp + 1].data() );
1980
1981 // Make some noise
1982 for( int nn = 0; nn < sz2Sided; ++nn )
1983 {
1984 N_n[nn] = normVar;
1985 N_nm[nn] = normVar;
1986 }
1987
1988 // Filter it
1989 // Set sqrt(PSD), just a pointer switch
1990 pfilt.psdSqrt( &sqrtNPSD[pp], tfreq[1] - tfreq[0] );
1991 nfilt.filter( N_n );
1992 nfilt.filter( N_nm );
1993
1994 // Normalize it
1995 realT Nactvar = 0.5 * ( math::vectorVariance( N_n ) + math::vectorVariance( N_nm ) );
1996 realT norm = sqrt( noiseVar[pp] / Nactvar );
1997 for( size_t q = 0; q < N_n.size(); ++q )
1998 N_n[q] *= norm;
1999 for( size_t q = 0; q < N_nm.size(); ++q )
2000 N_nm[q] *= norm;
2001
2002 // And move them to the Fourier domain
2003 fftF( Ntform1.data(), N_n.data() );
2004 fftF( Ntform2.data(), N_nm.data() );
2005
2006 // Apply the closed loop transfers
2007 for( size_t mm = 0; mm < tform1.size(); ++mm )
2008 {
2009 // Apply the augmented ETF to two time-series
2010 tform1lp[mm] = tform1[mm] * ETFlp[pp][mm] / scale;
2011 tform2lp[mm] = tform2[mm] * ETFlp[pp][mm] / scale;
2012
2013 Ntform1lp[mm] = Ntform1[mm] * NTFlp[pp][mm] / scale;
2014 Ntform2lp[mm] = Ntform2[mm] * NTFlp[pp][mm] / scale;
2015
2016 tform1[mm] *= ETFsi[pp][mm] / scale;
2017 tform2[mm] *= ETFsi[pp][mm] / scale;
2018
2019 Ntform1[mm] *= NTFsi[pp][mm] / scale;
2020 Ntform2[mm] *= NTFsi[pp][mm] / scale;
2021 }
2022
2023 // And make the speckle TS
2024 int m = fms[2 * pp].m;
2025 int n = fms[2 * pp].n;
2026
2027 //<<<<<<<<****** Transform back to the time domain.
2028 fftB( htsCorr[2 * pp].data(), tform1.data() );
2029 fftB( htsCorr[2 * pp + 1].data(), tform2.data() );
2030 fftB( N_n.data(), Ntform1.data() );
2031 fftB( N_nm.data(), Ntform2.data() );
2032
2033 for( int i = 0; i < hts[2 * pp].size(); ++i )
2034 {
2035 realT h1 = htsCorr[2 * pp][i] + N_n[i];
2036 realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2037
2038 spTS.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2039 spTS.image( i )( mnMax - m, mnMax - n ) = spTS.image( i )( mnMax + m, mnMax + n );
2040 }
2041
2042 fftB( htsCorr[2 * pp].data(), tform1lp.data() );
2043 fftB( htsCorr[2 * pp + 1].data(), tform2lp.data() );
2044 fftB( N_n.data(), Ntform1lp.data() );
2045 fftB( N_nm.data(), Ntform2lp.data() );
2046
2047 for( int i = 0; i < hts[2 * pp].size(); ++i )
2048 {
2049 realT h1 = htsCorr[2 * pp][i] + N_n[i];
2050 realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2051
2052 spTSlp.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2053 spTSlp.image( i )( mnMax - m, mnMax - n ) = spTSlp.image( i )( mnMax + m, mnMax + n );
2054 }
2055 }
2056
2057 /*********************************************************************/
2058 // 2.5) Calculate speckle PSD for each mode
2059 /*********************************************************************/
2060 std::vector<realT> speckAmp( spTS.planes() );
2061 std::vector<realT> speckAmplp( spTS.planes() );
2062
2063 for( size_t pp = 0; pp < nModes; ++pp )
2064 {
2065 int m = fms[2 * pp].m;
2066 int n = fms[2 * pp].n;
2067
2068 realT mn = 0;
2069 realT mnlp = 0;
2070 for( int i = 0; i < spTS.planes(); ++i )
2071 {
2072 speckAmp[i] = spTS.image( i )( mnMax + m, mnMax + n );
2073 speckAmplp[i] = spTSlp.image( i )( mnMax + m, mnMax + n );
2074
2075 mn += speckAmp[i];
2076 mnlp += speckAmplp[i];
2077 }
2078 mn /= speckAmp.size();
2079 mnlp /= speckAmplp.size();
2080
2081 // mean subtract
2082 for( int i = 0; i < speckAmp.size(); ++i )
2083 speckAmp[i] -= mn;
2084 for( int i = 0; i < speckAmplp.size(); ++i )
2085 speckAmplp[i] -= mnlp;
2086
2087 // Calculate PSD of the speckle amplitude
2088 avgPgram( tpgram, speckAmp );
2089 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2090 spPSDs[pp][nn] += tpgram[nn];
2091
2092 avgPgram( tpgram, speckAmplp );
2093 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2094 spPSDslp[pp][nn] += tpgram[nn];
2095 }
2096
2097 watcher.incrementAndOutputStatus();
2098 } // for(int zz=0; zz<lifetimeTrials; ++zz)
2099 } // #pragma omp parallel
2100
2101 std::vector<realT> spFreq( spPSDs[0].size() );
2102 for( size_t nn = 0; nn < spFreq.size(); ++nn )
2103 spFreq[nn] = tavgPgram[nn];
2104
2105 improc::eigenImage<realT> taus, tauslp;
2106 taus.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2107 tauslp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2108
2109 improc::eigenCube<realT> imc, imclp;
2110 std::vector<realT> clPSD;
2111
2112 if( writePSDs )
2113 {
2114 imc.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2115 imclp.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2116 clPSD.resize( sz2Sided );
2117 }
2118
2119 sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
2120 /*********************************************************************/
2121 // 3.0) Calculate lifetimes from the PSDs
2122 /*********************************************************************/
2123 for( size_t pp = 0; pp < nModes; ++pp )
2124 {
2125 spPSDs[pp][0] = spPSDs[pp][1]; // deal with under-estimated mean.
2126 spPSDslp[pp][0] = spPSDslp[pp][1]; // deal with under-estimated mean.
2127
2128 int m = fms[2 * pp].m;
2129 int n = fms[2 * pp].n;
2130
2131 realT var;
2132 if( writePSDs ) // Have to normalize the intensity for some reason if we want to use the PSDs
2133 {
2134 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2135 {
2136 spPSDs[pp][nn] /= lifetimeTrials;
2137 }
2138
2139 for( size_t nn = 0; nn < sz2Sided; ++nn )
2140 {
2141 clPSD[nn] =
2142 opdPSD[pp][nn] * norm( ETFsi[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFsi[pp][nn] );
2143 }
2144
2145 var = sigproc::psdVar( freq2sided, clPSD );
2146
2147 realT pvar = sigproc::psdVar( spFreq, spPSDs[pp] );
2148
2149 for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2150 {
2151 spPSDs[pp][nn] *= var / pvar;
2152 imc.image( nn )( mnMax + m, mnMax + n ) = spPSDs[pp][nn];
2153 imc.image( nn )( mnMax - m, mnMax - n ) = spPSDs[pp][nn];
2154 }
2155
2156 // lp
2157 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2158 {
2159 spPSDslp[pp][nn] /= lifetimeTrials;
2160 }
2161
2162 for( size_t nn = 0; nn < sz2Sided; ++nn )
2163 {
2164 clPSD[nn] =
2165 opdPSD[pp][nn] * norm( ETFlp[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFlp[pp][nn] );
2166 }
2167
2168 var = sigproc::psdVar( freq2sided, clPSD );
2169
2170 pvar = sigproc::psdVar( spFreq, spPSDslp[pp] );
2171
2172 for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2173 {
2174 spPSDslp[pp][nn] *= var / pvar;
2175 imclp.image( nn )( mnMax + m, mnMax + n ) = spPSDslp[pp][nn];
2176 imclp.image( nn )( mnMax - m, mnMax - n ) = spPSDslp[pp][nn];
2177 }
2178 }
2179
2180 var = sigproc::psdVar( spFreq, spPSDs[pp] );
2181
2182 realT T = ( 1.0 / ( spFreq[1] - spFreq[0] ) ) * 10;
2183 realT error;
2184 realT tau = pvm( error, spFreq, spPSDs[pp], T ) * ( T ) / var;
2185 taus( mnMax + m, mnMax + n ) = tau;
2186 taus( mnMax - m, mnMax - n ) = tau;
2187
2188 var = sigproc::psdVar( spFreq, spPSDslp[pp] );
2189
2190 tau = pvm( error, spFreq, spPSDslp[pp], T ) * ( T ) / var;
2191 tauslp( mnMax + m, mnMax + n ) = tau;
2192 tauslp( mnMax - m, mnMax - n ) = tau;
2193 }
2194
2195 /*********************************************************************/
2196 // 4.0) Write the results to disk
2197 /*********************************************************************/
2198 fn = std::format("{}/speckleLifetimes_{}_si.fits",dir, mags[s]);
2199 //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2200 ff.write( fn, taus );
2201
2202 fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]);
2203 //dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2204 ff.write( fn, tauslp );
2205
2206 if( writePSDs )
2207 {
2208 fn = std::format("{}/specklePSDs_{}_si.fits",dir, mags[s]);
2209 //dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2210 ff.write( fn, imc );
2211
2212 fn = std::format("{}/speckleLifetimes_{}_lp.fits",dir, mags[s]);
2213 //dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2214 ff.write( fn, imclp );
2215 }
2216
2217 } // mags
2218
2219 return 0;
2220}
2221
2222template <typename realT, typename aosysT>
2223int fourierTemporalPSD<realT, aosysT>::getGridFreq( std::vector<realT> &freq, const std::string &dir )
2224{
2225 std::string fn;
2226 fn = dir + "/psds/freq.binv";
2227 return ioutils::readBinVector( freq, fn );
2228}
2229
2230template <typename realT, typename aosysT>
2231int fourierTemporalPSD<realT, aosysT>::getGridPSD( std::vector<realT> &psd, const std::string &dir, int m, int n )
2232{
2233 std::string fn;
2234 fn = std::format("{}/psds/psd_{}_{}.binv",dir,m,n);
2235 //dir + "/psds/psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2236 return ioutils::readBinVector( psd, fn );
2237}
2238
2239template <typename realT, typename aosysT>
2240int fourierTemporalPSD<realT, aosysT>::getGridPSD(
2241 std::vector<realT> &freq, std::vector<realT> &psd, const std::string &dir, int m, int n )
2242{
2243 int rv = getGridFreq( freq, dir );
2244 if( rv < 0 )
2245 return rv;
2246
2247 return getGridPSD( psd, dir, m, n );
2248}
2249
2250/// Worker function for GSL Integration for the basic sin/cos Fourier modes.
2251/** \ingroup mxAOAnalytic
2252 */
2253template <typename realT, typename aosysT>
2254realT F_basic( realT kv, void *params )
2255{
2256 fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2257
2258 realT f = Fp->m_f;
2259 realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2260
2261 realT D = Fp->m_aosys->D();
2262 realT m = Fp->m_m;
2263 realT n = Fp->m_n;
2264 int p = Fp->m_p;
2265
2266 realT ku = f / v_wind;
2267
2268 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2269 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2270
2271 realT Q1 = math::func::jinc( math::pi<realT>() * D * kp );
2272
2273 realT Q2 = math::func::jinc( math::pi<realT>() * D * kpp );
2274
2275 realT Q = ( Q1 + p * Q2 );
2276
2277 realT P =
2278 Fp->m_aosys->psd( Fp->m_aosys->atm, Fp->_layer_i, sqrt( pow( ku, 2 ) + pow( kv, 2 ) ), Fp->m_aosys->secZeta() );
2279
2280 return P * Q * Q;
2281}
2282
2283template <typename realT>
2284void turbBoilCubic(
2285 realT &a, realT &b, realT &c, realT &d, const realT &kv, const realT &f, const realT &Vu, const realT &f0, int pm )
2286{
2287 a = Vu * Vu * Vu;
2288 b = -( 3 * Vu * Vu * f + pm * f0 * f0 * f0 );
2289 c = 3 * f * f * Vu;
2290 d = -( f * f * f + pm * f0 * f0 * f0 * kv * kv );
2291}
2292
2293/// Worker function for GSL Integration for the modified Fourier modes.
2294/** \ingroup mxAOAnalytic
2295 */
2296template <typename realT, typename aosysT>
2297realT F_mod( realT kv, void *params )
2298{
2299 fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2300
2301 realT f = Fp->m_f;
2302 realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2303
2304 realT D = Fp->m_aosys->D();
2305 realT m = Fp->m_m;
2306 realT n = Fp->m_n;
2307
2308 realT f0 = Fp->m_f0;
2309
2310 realT ku;
2311 if( f0 == 0 )
2312 {
2313 ku = f / v_wind;
2314
2315 if( Fp->m_spatialFilter )
2316 {
2317 // de-rotate the spatial frequency vector back to pupil coordinates
2318 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2319 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2320 // Return if spatially filtered
2321 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2322 return 0;
2323
2324 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2325 return 0;
2326 }
2327
2328 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2329 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2330
2331 realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
2332
2333 realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2334
2335 realT QQ = 2 * ( Jp * Jp + Jm * Jm );
2336
2337 realT P = Fp->m_aosys->psd( Fp->m_aosys->atm,
2338 Fp->_layer_i,
2339 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2340 Fp->m_aosys->lam_sci(),
2341 Fp->m_aosys->lam_wfs(),
2342 Fp->m_aosys->secZeta() );
2343
2344 return P * QQ;
2345 }
2346 else
2347 {
2348 realT a, b, c, d, p, q;
2349
2350 turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, 1 );
2351 mx::math::cubicDepressed( p, q, a, b, c, d );
2352 realT t = mx::math::cubicRealRoot( p, q );
2353
2354 ku = t - b / ( 3 * a );
2355
2356 if( Fp->m_spatialFilter )
2357 {
2358 // de-rotate the spatial frequency vector back to pupil coordinates
2359 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2360 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2361 // Return if spatially filtered
2362 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2363 return 0;
2364
2365 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2366 return 0;
2367 }
2368
2369 realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2370 realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2371
2372 realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
2373
2374 realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2375
2376 realT QQ = 2 * ( Jp * Jp + Jm * Jm );
2377
2378 realT P1 = Fp->m_aosys->psd( Fp->m_aosys->atm,
2379 Fp->_layer_i,
2380 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2381 Fp->m_aosys->lam_sci(),
2382 Fp->m_aosys->lam_wfs(),
2383 Fp->m_aosys->secZeta() );
2384
2385 P1 *= QQ;
2386
2387 turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, -1 );
2388 mx::math::cubicDepressed( p, q, a, b, c, d );
2389 t = mx::math::cubicRealRoot( p, q );
2390
2391 ku = t - b / ( 3 * a );
2392
2393 if( Fp->m_spatialFilter )
2394 {
2395 // de-rotate the spatial frequency vector back to pupil coordinates
2396 realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2397 realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2398 // Return if spatially filtered
2399 if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2400 return 0;
2401
2402 if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2403 return 0;
2404 }
2405
2406 kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2407 kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2408
2409 Jp = math::func::jinc( math::pi<realT>() * D * kp );
2410
2411 Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2412
2413 QQ = 2 * ( Jp * Jp + Jm * Jm );
2414
2415 realT P2 = Fp->m_aosys->psd( Fp->m_aosys->atm,
2416 Fp->_layer_i,
2417 sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2418 Fp->m_aosys->lam_sci(),
2419 Fp->m_aosys->lam_wfs(),
2420 Fp->m_aosys->secZeta() );
2421
2422 P2 *= QQ;
2423
2424 return 0.5*(P1 + P2);
2425 }
2426}
2427
2428
2429/*extern template
2430struct fourierTemporalPSD<float, aoSystem<float, vonKarmanSpectrum<float>, std::ostream>>;*/
2431
2432extern template struct fourierTemporalPSD<double, aoSystem<double, vonKarmanSpectrum<double>, std::ostream>>;
2433
2434/*
2435extern template
2436struct fourierTemporalPSD<long double, aoSystem<long double, vonKarmanSpectrum<long double>, std::ostream>>;
2437
2438#ifdef HASQUAD
2439extern template
2440struct fourierTemporalPSD<__float128, aoSystem<__float128, vonKarmanSpectrum<__float128>, std::ostream>>;
2441#endif
2442*/
2443
2444} // namespace analysis
2445} // namespace AO
2446} // namespace mx
2447
2448#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:53
mx::fits::fitsFile::read
error_t read(dataT *data)
Read the contents of the FITS file into an array.
Definition fitsFile.hpp:992
mx::fits::fitsFile::write
error_t write(const dataT *im, int d1, int d2, int d3, fitsHeader< verboseT > *head)
Write the contents of a raw array to the FITS file.
Definition fitsFile.hpp:1587
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:240
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::invalidarg
@ invalidarg
An argument was invalid.
mx::error_t::error
@ error
A general error has occurred.
mx::internal::mxlib_error_report
error_t mxlib_error_report(const error_t &code, const std::string &expl, const std::source_location &loc=std::source_location::current())
Print a report to stderr given an mxlib error_t code and explanation and return the code.
Definition error.hpp:331
mx::ioutils::createDirectories
error_t createDirectories(const std::string &path)
Create a directory or directories.
Definition fileUtils.cpp:54
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:542
mx::sigproc::makeFourierModeFreqs_Rect
int makeFourierModeFreqs_Rect(std::vector< fourierModeDef > &spf, int N)
Generate a rectangular spatial frequency grid.
Definition fourierModes.hpp:593
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:2254
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:2297
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:834
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:893
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:138
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:453
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::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 mxlib.hpp:37
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:893
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:1143
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:962
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:851
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:931
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)
*‍/
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:2223
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:2231
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:1531
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.