Line data Source code
1 : /** \file fourierTemporalPSD.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Calculation of the temporal PSD of Fourier modes.
4 : * \ingroup mxAO_analysis_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 <algorithm>
31 : #include <atomic>
32 : #include <cmath>
33 : #include <iostream>
34 : #include <fstream>
35 : #include <limits>
36 : #include <map>
37 : #include <memory>
38 : #include <mutex>
39 :
40 : #include <sys/stat.h>
41 :
42 : #include <gsl/gsl_integration.h>
43 : #include <gsl/gsl_errno.h>
44 :
45 : #include <Eigen/Dense>
46 :
47 : #include "../../mxlib.hpp"
48 : #include "../../math/constants.hpp"
49 : #include "../../math/floatUtils.hpp"
50 : #include "../../math/func/jinc.hpp"
51 : #include "../../math/func/airyPattern.hpp"
52 : #include "../../math/vectorUtils.hpp"
53 : #include "../../ioutils/fits/fitsFile.hpp"
54 : #include "../../sigproc/fourierModes.hpp"
55 : #include "../../sigproc/psdVarMean.hpp"
56 : #include "../../ioutils/stringUtils.hpp"
57 : #include "../../ioutils/readColumns.hpp"
58 : #include "../../ioutils/binVector.hpp"
59 : #include "../../ioutils/fileUtils.hpp"
60 :
61 : #include "../../ipc/ompLoopWatcher.hpp"
62 :
63 : #include "aoSystem.hpp"
64 : #include "aoPSDs.hpp"
65 : #include "wfsNoisePSD.hpp"
66 : #include "clAOLinearPredictor.hpp"
67 : #include "clGainOpt.hpp"
68 : #include "varmapToImage.hpp"
69 : #include "speckleAmpPSD.hpp"
70 :
71 : #include "aoConstants.hpp"
72 :
73 : namespace mx
74 : {
75 : namespace AO
76 : {
77 : namespace analysis
78 : {
79 :
80 : #ifndef WSZ
81 :
82 : /** \def WSZ
83 : * \brief Size of the GSL integration workspace
84 : */
85 : #define WSZ 100000
86 :
87 : #endif
88 :
89 : enum basis : unsigned int
90 : {
91 : basic, ///< The basic sine and cosine Fourier modes
92 : modified ///< The modified Fourier basis from \cite males_guyon_2018
93 : };
94 :
95 : /// Policy for handling GSL quadrature non-convergence statuses.
96 : enum class fourierTemporalPSDPolicy
97 : {
98 : permissive, ///< Retain the best finite approximation and record the status.
99 : strict ///< Record every status and return an error if any integration does not converge.
100 : };
101 :
102 : /// Aggregated GSL quadrature diagnostics for a Fourier temporal PSD calculation.
103 : template <typename realT>
104 : struct fourierTemporalPSDReport
105 : {
106 : /// Summary of one GSL status code.
107 : struct statusSummary
108 : {
109 : size_t count{ 0 }; ///< Number of occurrences of this status.
110 : std::map<size_t, size_t> countByLayer; ///< Number of occurrences in each atmospheric layer.
111 : realT maximumAbsoluteError{ 0 }; ///< Largest GSL absolute-error estimate.
112 : realT maximumToleranceRatio{ 0 }; ///< Largest error estimate relative to the requested tolerance.
113 : size_t worstLayer{ 0 }; ///< Layer containing the largest tolerance ratio.
114 : realT worstFrequency{ 0 }; ///< Frequency containing the largest tolerance ratio.
115 : };
116 :
117 : size_t integrationsAttempted{ 0 }; ///< Total number of quadrature calls.
118 : size_t integrationsConverged{ 0 }; ///< Number of quadrature calls returning `GSL_SUCCESS`.
119 : std::map<int, statusSummary> gslStatus; ///< Summaries keyed by the raw GSL status code.
120 :
121 : /// Reset all accumulated diagnostics.
122 : void clear();
123 :
124 : /// Record one quadrature result.
125 : void record( int status, /**< [in] raw GSL status code */
126 : size_t layer, /**< [in] atmospheric layer index */
127 : realT frequency, /**< [in] temporal frequency */
128 : realT result, /**< [in] quadrature result */
129 : realT absoluteError, /**< [in] GSL absolute-error estimate */
130 : realT absoluteTolerance, /**< [in] requested absolute tolerance */
131 : realT relativeTolerance /**< [in] requested relative tolerance */ );
132 :
133 : /// Merge another report into this report.
134 : void merge( const fourierTemporalPSDReport &other /**< [in] report to merge */ );
135 :
136 : /// Return the total number of non-successful integrations.
137 : [[nodiscard]] size_t failureCount() const;
138 :
139 : /// Write a human-readable summary of the accumulated quadrature diagnostics.
140 : void write( std::ostream &output /**< [out] stream receiving the summary */ ) const;
141 : };
142 :
143 : /// \cond fourierTemporalPSD_detail
144 : namespace fourierTemporalPSD_detail
145 : {
146 :
147 : /// Function type used to allocate a GSL integration workspace.
148 : using gslWorkspaceAllocator = gsl_integration_workspace *(*)( size_t );
149 :
150 : /// Deleter providing RAII ownership for a GSL integration workspace.
151 : struct gslWorkspaceDeleter
152 : {
153 : /// Free an allocated GSL integration workspace.
154 7 : void operator()( gsl_integration_workspace *workspace /**< [in] workspace to free */ ) const noexcept
155 : {
156 7 : if( workspace != nullptr )
157 : {
158 7 : gsl_integration_workspace_free( workspace );
159 : }
160 7 : }
161 : };
162 :
163 : /// Unique ownership handle for a GSL integration workspace.
164 : using gslWorkspacePtr = std::unique_ptr<gsl_integration_workspace, gslWorkspaceDeleter>;
165 :
166 : /// Return whether a GSL status represents a potentially usable non-converged approximation.
167 12 : inline bool isConvergenceStatus( int status )
168 : {
169 12 : return status == GSL_EMAXITER || status == GSL_EROUND || status == GSL_ESING || status == GSL_EDIVERGE;
170 : }
171 :
172 : /// Convert a fatal GSL status to an mxlib status.
173 4 : inline error_t gslStatusToError( int status )
174 : {
175 4 : if( status == GSL_ENOMEM )
176 : {
177 1 : return error_t::allocerr;
178 : }
179 :
180 3 : if( status == GSL_EDOM || status == GSL_EINVAL )
181 : {
182 2 : return error_t::invalidconfig;
183 : }
184 :
185 1 : return error_t::liberr;
186 : }
187 :
188 : /// Apply the configured non-convergence policy to one GSL status.
189 116 : inline error_t applyPolicy( int status, fourierTemporalPSDPolicy policy )
190 : {
191 116 : if( status == GSL_SUCCESS )
192 : {
193 104 : return error_t::noerror;
194 : }
195 :
196 12 : if( isConvergenceStatus( status ) )
197 : {
198 8 : return policy == fourierTemporalPSDPolicy::permissive ? error_t::noerror : error_t::liberr;
199 : }
200 :
201 4 : return gslStatusToError( status );
202 : }
203 :
204 : /// Mutex serializing scoped changes to GSL's process-global error handler.
205 12 : inline std::mutex &gslErrorHandlerMutex()
206 : {
207 : static std::mutex mutex;
208 12 : return mutex;
209 : }
210 :
211 : /// Disable the GSL error handler for a complete top-level PSD calculation and restore it on exit.
212 : /** The mutex serializes handler changes made by this implementation. Unrelated code cannot be protected from GSL's
213 : * process-global handler state unless it coordinates with the same mutex.
214 : */
215 : class scopedGslErrorHandlerOff
216 : {
217 : public:
218 : /// Lock handler management and retain the previously installed handler.
219 12 : scopedGslErrorHandlerOff() : m_lock( gslErrorHandlerMutex() ), m_previous( gsl_set_error_handler_off() )
220 : {
221 12 : }
222 :
223 : /// Disallow copying ownership of the saved handler.
224 : scopedGslErrorHandlerOff( const scopedGslErrorHandlerOff & ) = delete;
225 :
226 : /// Disallow copy assignment of the handler guard.
227 : scopedGslErrorHandlerOff &operator=( const scopedGslErrorHandlerOff & ) = delete;
228 :
229 : /// Restore the previously installed handler before releasing the lock.
230 12 : ~scopedGslErrorHandlerOff()
231 : {
232 12 : static_cast<void>( gsl_set_error_handler( m_previous ) );
233 12 : }
234 :
235 : private:
236 : std::unique_lock<std::mutex> m_lock; ///< Lock held while the handler is disabled.
237 : gsl_error_handler_t *m_previous{ nullptr }; ///< Handler restored on destruction.
238 : };
239 :
240 : } // namespace fourierTemporalPSD_detail
241 : /// \endcond
242 :
243 : template <typename realT>
244 29 : void fourierTemporalPSDReport<realT>::clear()
245 : {
246 29 : integrationsAttempted = 0;
247 29 : integrationsConverged = 0;
248 29 : gslStatus.clear();
249 29 : }
250 :
251 : template <typename realT>
252 122 : void fourierTemporalPSDReport<realT>::record( int status,
253 : size_t layer,
254 : realT frequency,
255 : realT result,
256 : realT absoluteError,
257 : realT absoluteTolerance,
258 : realT relativeTolerance )
259 : {
260 122 : ++integrationsAttempted;
261 122 : if( status == GSL_SUCCESS )
262 : {
263 119 : ++integrationsConverged;
264 119 : return;
265 : }
266 :
267 3 : statusSummary &summary = gslStatus[status];
268 3 : ++summary.count;
269 3 : ++summary.countByLayer[layer];
270 3 : summary.maximumAbsoluteError = std::max( summary.maximumAbsoluteError, std::abs( absoluteError ) );
271 :
272 3 : const realT requestedTolerance = std::max( std::abs( absoluteTolerance ), std::abs( relativeTolerance * result ) );
273 3 : const realT toleranceRatio = requestedTolerance > 0 ? std::abs( absoluteError ) / requestedTolerance
274 0 : : std::numeric_limits<realT>::infinity();
275 3 : if( toleranceRatio >= summary.maximumToleranceRatio )
276 : {
277 3 : summary.maximumToleranceRatio = toleranceRatio;
278 3 : summary.worstLayer = layer;
279 3 : summary.worstFrequency = frequency;
280 : }
281 : }
282 :
283 : template <typename realT>
284 5 : void fourierTemporalPSDReport<realT>::merge( const fourierTemporalPSDReport &other )
285 : {
286 5 : integrationsAttempted += other.integrationsAttempted;
287 5 : integrationsConverged += other.integrationsConverged;
288 :
289 6 : for( const auto &[status, otherSummary] : other.gslStatus )
290 : {
291 1 : statusSummary &summary = gslStatus[status];
292 1 : summary.count += otherSummary.count;
293 2 : for( const auto &[layer, count] : otherSummary.countByLayer )
294 : {
295 1 : summary.countByLayer[layer] += count;
296 : }
297 1 : summary.maximumAbsoluteError = std::max( summary.maximumAbsoluteError, otherSummary.maximumAbsoluteError );
298 1 : if( otherSummary.maximumToleranceRatio >= summary.maximumToleranceRatio )
299 : {
300 1 : summary.maximumToleranceRatio = otherSummary.maximumToleranceRatio;
301 1 : summary.worstLayer = otherSummary.worstLayer;
302 1 : summary.worstFrequency = otherSummary.worstFrequency;
303 : }
304 : }
305 5 : }
306 :
307 : template <typename realT>
308 1 : size_t fourierTemporalPSDReport<realT>::failureCount() const
309 : {
310 1 : return integrationsAttempted - integrationsConverged;
311 : }
312 :
313 : template <typename realT>
314 1 : void fourierTemporalPSDReport<realT>::write( std::ostream &output ) const
315 : {
316 1 : output << "GSL quadrature: " << integrationsConverged << '/' << integrationsAttempted << " converged\n";
317 3 : for( const auto &[status, summary] : gslStatus )
318 : {
319 2 : output << " " << gsl_strerror( status ) << " (" << status << "): " << summary.count << ", max absolute error "
320 2 : << summary.maximumAbsoluteError << ", max tolerance ratio " << summary.maximumToleranceRatio
321 2 : << " at layer " << summary.worstLayer << ", frequency " << summary.worstFrequency << ", layers {";
322 2 : bool firstLayer = true;
323 5 : for( const auto &[layer, count] : summary.countByLayer )
324 : {
325 3 : if( !firstLayer )
326 : {
327 1 : output << ", ";
328 : }
329 3 : output << layer << ": " << count;
330 3 : firstLayer = false;
331 : }
332 2 : output << "}\n";
333 : }
334 1 : }
335 :
336 : // Forward declaration
337 : template <typename realT, typename aosysT>
338 : realT F_basic( realT kv, void *params );
339 :
340 : // Forward declaration
341 : template <typename realT, typename aosysT>
342 : realT F_mod( realT kv, void *params );
343 :
344 : /// Class to manage the calculation of temporal PSDs of the Fourier modes in atmospheric turbulence.
345 : /** Works with both basic (sines/cosines) and modified Fourier modes.
346 : *
347 : * \tparam realT is a real floating point type for calculations. Currently must be double due to gsl_integration.
348 : * \tparam aosysT is an AO system type, usually of type ao_system.
349 : *
350 : * \todo Split off the integration parameters in a separate structure.
351 : * \todo once integration parameters are in a separate structure, make this a class with protected members.
352 : * \ingroup mxAOAnalytic
353 : */
354 : template <typename _realT, typename aosysT>
355 : struct fourierTemporalPSD
356 : {
357 : /// The type for arithmetic
358 : typedef _realT realT;
359 :
360 : /// The complex type for arithmetic
361 : typedef std::complex<realT> complexT;
362 :
363 : /// Quadrature report type used by this specialization.
364 : typedef fourierTemporalPSDReport<realT> reportT;
365 :
366 : /// Pointer to an AO system structure.
367 : aosysT *m_aosys{ nullptr };
368 :
369 : realT m_f{ 0 }; ///< the current temporal frequency
370 : realT m_m{ 0 }; ///< the spatial frequency m index
371 : realT m_n{ 0 }; ///< the spatial frequency n index
372 : realT m_cq{ 0 }; ///< The cosine of the wind direction
373 : realT m_sq{ 0 }; ///< The sine of the wind direction
374 : realT m_spatialFilter{ false }; ///< Flag indicating if a spatial filter is applied
375 :
376 : bool m_strehlOG{ false };
377 : bool m_uncorrectedOG{ false };
378 :
379 : realT m_f0{ 0 }; ///< the Berdja boiling parameter
380 :
381 : int m_p{ 1 }; ///< The parity of the mode, +/- 1. If _useBasis==MXAO_FTPSD_BASIS_BASIC then +1 indicates cosine, -1
382 : ///< indicates sine.
383 : int _layer_i; ///< The index of the current layer.
384 :
385 : int _useBasis; ///< Set to MXAO_FTPSD_BASIS_BASIC/MODIFIED/PROJECTED_* to use the basic sin/cos modes, the modified
386 : ///< Fourier modes, or a projection of them.
387 :
388 : protected:
389 : /// Unique ownership of the GSL integration workspace used by worker instances.
390 : fourierTemporalPSD_detail::gslWorkspacePtr m_workspace;
391 :
392 : /// Allocation function used when a worker lazily creates its GSL workspace.
393 : fourierTemporalPSD_detail::gslWorkspaceAllocator m_workspaceAllocator{ gsl_integration_workspace_alloc };
394 :
395 : public:
396 : realT _absTol; ///< The absolute tolerance to use in the GSL integrator
397 : realT _relTol; ///< The relative tolerance to use in the GSL integrator
398 :
399 : int m_mode_i; ///< Projected basis mode index
400 :
401 : Eigen::Array<realT, -1, -1> m_modeCoeffs; ///< Coeeficients of the projection onto the Fourier modes
402 : realT m_minCoeffVal;
403 :
404 : std::vector<realT> Jps;
405 : std::vector<realT> Jms;
406 : std::vector<int> ps;
407 : std::vector<realT> ms;
408 : std::vector<realT> ns;
409 :
410 0 : void initProjection()
411 : {
412 0 : Jps.resize( m_modeCoeffs.cols() );
413 0 : Jms.resize( m_modeCoeffs.cols() );
414 0 : ps.resize( m_modeCoeffs.cols() );
415 0 : ms.resize( m_modeCoeffs.cols() );
416 0 : ns.resize( m_modeCoeffs.cols() );
417 :
418 0 : for( int i = 0; i < m_modeCoeffs.cols(); ++i )
419 : {
420 : int m, n, p;
421 0 : sigproc::fourierModeCoordinates( m, n, p, i );
422 0 : ps[i] = p;
423 0 : ms[i] = m;
424 0 : ns[i] = n;
425 : }
426 0 : }
427 :
428 : public:
429 : /// Default c'tor
430 : fourierTemporalPSD();
431 :
432 : /// Disallow copying unique workspace ownership.
433 : fourierTemporalPSD( const fourierTemporalPSD & ) = delete;
434 :
435 : /// Disallow copy assignment of unique workspace ownership.
436 : fourierTemporalPSD &operator=( const fourierTemporalPSD & ) = delete;
437 :
438 : /// Move workspace ownership and evaluator state.
439 1 : fourierTemporalPSD( fourierTemporalPSD && ) noexcept = default;
440 :
441 : /// Move-assign workspace ownership and evaluator state.
442 : fourierTemporalPSD &operator=( fourierTemporalPSD && ) noexcept = default;
443 :
444 : /// Release owned resources.
445 48 : ~fourierTemporalPSD() = default;
446 :
447 : protected:
448 : /// Construct with a custom workspace allocator.
449 : explicit fourierTemporalPSD(
450 : fourierTemporalPSD_detail::gslWorkspaceAllocator allocator /**< [in] workspace allocation function */ );
451 :
452 : /// Initialize parameters to default values.
453 : void initialize();
454 :
455 : /// Allocate the worker workspace if it is not already available.
456 : error_t allocateWorkspace();
457 :
458 : /// Validate state and arguments shared by single- and multilayer calculations.
459 : error_t validatePsdInputs( const std::vector<realT> &PSD, /**< [in] output storage to validate */
460 : const std::vector<realT> &freq, /**< [in] temporal-frequency grid */
461 : realT m, /**< [in] first spatial-frequency index */
462 : realT n, /**< [in] second spatial-frequency index */
463 : int p, /**< [in] Fourier-mode parity */
464 : realT fmax, /**< [in] maximum exactly integrated frequency */
465 : int layer_i, /**< [in] layer index, or -1 to validate all layers */
466 : fourierTemporalPSDPolicy policy /**< [in] non-convergence policy */ );
467 :
468 : /// Validate the configured atmosphere and optionally a requested layer.
469 : error_t validateAtmosphere( int layer_i /**< [in] layer index, or -1 to validate all layers */ );
470 :
471 : public:
472 : /** \name GSL Integration Tolerances
473 : * For good results it seems that absolute tolerance (absTol) needs to be 1e-10. Lower tolerances cause some
474 : * frequencies to drop out, etc. Relative tolerance (relTol) seems to be less sensitive, and 1e-4 works on cases
475 : * tested as of 1 Jan, 2017.
476 : *
477 : * See the documentation for the GSL Library integrators at
478 : * (https://www.gnu.org/software/gsl/manual/htmlm_node/QAGI-adaptive-integration-on-infinite-intervals.html)
479 : * @{
480 : */
481 :
482 : /// Set absolute tolerance
483 : /**
484 : * \param at is the new absolute tolerance.
485 : */
486 : void absTol( realT at );
487 :
488 : /// Get the current absolute tolerance
489 : /**
490 : * \returns _absTol
491 : */
492 : realT absTol();
493 :
494 : /// Set relative tolerance
495 : /**
496 : * \param rt is the new relative tolerance.
497 : */
498 : void relTol( realT rt );
499 :
500 : /// Get the current relative tolerance
501 : /**
502 : * \returns _relTol
503 : */
504 : realT relTol();
505 :
506 : ///@}
507 :
508 : /// Determine the frequency of the highest V-dot-k peak
509 : /**
510 : * \param m the spatial frequency u index
511 : * \param n the spatial frequency v index
512 : *
513 : * \return the frequency of the fastest peak
514 : */
515 : realT fastestPeak( int m, int n );
516 :
517 : protected:
518 : /// Calculate a single-layer temporal PSD while the caller manages the GSL error handler.
519 : error_t singleLayerPSDImpl( std::vector<realT> &PSD, /**< [out] calculated PSD */
520 : std::vector<realT> &freq, /**< [in] temporal-frequency grid */
521 : realT m, /**< [in] first spatial-frequency index */
522 : realT n, /**< [in] second spatial-frequency index */
523 : int layer_i, /**< [in] atmospheric-layer index */
524 : int p, /**< [in] Fourier-mode parity */
525 : realT fmax, /**< [in] maximum exactly integrated frequency */
526 : reportT &report, /**< [out] accumulated quadrature report */
527 : fourierTemporalPSDPolicy policy /**< [in] non-convergence policy */ );
528 :
529 : public:
530 : /// Calculate the temporal PSD for a Fourier mode for a single layer.
531 : /** `PSD` and `freq` must have the same nonzero size. Frequencies must be finite, nonnegative, and strictly
532 : * increasing. The AO system, integration controls, requested layer, and atmosphere are validated before
533 : * calculation. A precondition or allocation failure leaves `PSD` unchanged and clears `report` when supplied.
534 : *
535 : * When extending beyond `fmax`, up to the last 50 exactly integrated bins are averaged after projection to the
536 : * first tail frequency. If fewer than 50 exact bins are available, all available exact bins are used. At least one
537 : * exact bin is required to initialize the tail.
538 : *
539 : * In permissive mode, finite best approximations returned with `GSL_EMAXITER`, `GSL_EROUND`, `GSL_ESING`, or
540 : * `GSL_EDIVERGE` are retained and summarized in `report`. In strict mode the calculation continues to characterize
541 : * all such failures but returns `error_t::liberr` and the output PSD must be discarded.
542 : *
543 : * \returns `error_t::noerror` on success, an argument/configuration error for invalid inputs, or
544 : * `error_t::liberr` when strict quadrature handling detects non-convergence.
545 : */
546 : error_t singleLayerPSD(
547 : std::vector<realT> &PSD, /**< [out] calculated PSD */
548 : std::vector<realT> &freq, /**< [in] temporal-frequency grid */
549 : realT m, /**< [in] first spatial-frequency index */
550 : realT n, /**< [in] second spatial-frequency index */
551 : int layer_i, /**< [in] atmospheric-layer index */
552 : int p, /**< [in] Fourier-mode parity */
553 : realT fmax = 0, /**< [in] maximum exactly integrated frequency, or 0 for the grid maximum */
554 : reportT *report = nullptr, /**< [out] optional quadrature report */
555 : fourierTemporalPSDPolicy policy = fourierTemporalPSDPolicy::permissive /**< [in] non-convergence policy */ );
556 :
557 : ///\cond multilayerm_parallel
558 : // Conditional to exclude from Doxygen.
559 :
560 : protected:
561 : // Type to allow overloading of the multiLayerPSD workers based on whether they are parallelized or not.
562 : template <bool m_parallel>
563 : struct isParallel
564 : {
565 : };
566 :
567 : // Parallelized version of multiLayerPSD, with OMP directives.
568 : error_t multiLayerPSDImpl( std::vector<realT> &PSD, /**< [out] calculated PSD */
569 : std::vector<realT> &freq, /**< [in] temporal-frequency grid */
570 : realT m, /**< [in] first spatial-frequency index */
571 : realT n, /**< [in] second spatial-frequency index */
572 : int p, /**< [in] Fourier-mode parity */
573 : realT fmax, /**< [in] maximum exactly integrated frequency */
574 : reportT &report, /**< [out] accumulated quadrature report */
575 : fourierTemporalPSDPolicy policy, /**< [in] non-convergence policy */
576 : isParallel<true> parallel /**< [in] parallel dispatch tag */ );
577 :
578 : // Non-Parallelized version of multiLayerPSD, without OMP directives.
579 : error_t multiLayerPSDImpl( std::vector<realT> &PSD, /**< [out] calculated PSD */
580 : std::vector<realT> &freq, /**< [in] temporal-frequency grid */
581 : realT m, /**< [in] first spatial-frequency index */
582 : realT n, /**< [in] second spatial-frequency index */
583 : int p, /**< [in] Fourier-mode parity */
584 : realT fmax, /**< [in] maximum exactly integrated frequency */
585 : reportT &report, /**< [out] accumulated quadrature report */
586 : fourierTemporalPSDPolicy policy, /**< [in] non-convergence policy */
587 : isParallel<false> parallel /**< [in] sequential dispatch tag */ );
588 :
589 : ///\endcond
590 :
591 : public:
592 : /// Calculate the temporal PSD for a Fourier mode in a multi-layer model.
593 : /** `PSD` and `freq` must have the same nonzero size. Frequencies must be finite, nonnegative, and strictly
594 : * increasing. The AO system, integration controls, and complete atmosphere are validated before calculation. A
595 : * precondition failure leaves `PSD` unchanged and clears `report` when supplied.
596 : *
597 : * \tparam parallel controls whether layers are calculated in parallel. Default is true. Set to false if this is
598 : * called inside a parallelized loop, as in \ref makePSDGrid.
599 : *
600 : * In permissive mode, recognized convergence failures are retained and summarized in `report`. Strict mode returns
601 : * `error_t::liberr` if any layer has such a failure; the output PSD is incomplete and must be discarded whenever
602 : * this function returns an error.
603 : *
604 : * \returns `error_t::noerror` on success, or the first layer error in atmospheric-layer order.
605 : */
606 : template <bool parallel = true>
607 : error_t multiLayerPSD(
608 : std::vector<realT> &PSD, /**< [out] calculated PSD */
609 : std::vector<realT> &freq, /**< [in] temporal-frequency grid */
610 : realT m, /**< [in] first spatial-frequency index */
611 : realT n, /**< [in] second spatial-frequency index */
612 : int p, /**< [in] Fourier-mode parity */
613 : realT fmax = 0, /**< [in] maximum exactly integrated frequency, or 0 for the default cutoff */
614 : reportT *report = nullptr, /**< [out] optional quadrature report */
615 : fourierTemporalPSDPolicy policy = fourierTemporalPSDPolicy::permissive /**< [in] non-convergence policy */ );
616 :
617 : /// Calculate PSDs over a grid of spatial frequencies.
618 : /** The grid of spatial frequencies is square, set by the maximum value of m and n.
619 : *
620 : * The PSDs are written as mx::binVector binary files to a directory. We do not use FITS since
621 : * this adds overhead and cfitisio handles parallelization poorly due to the limitation on number of created file
622 : * pointers.
623 : *
624 : * Inputs and AO-system state are validated before any output is created. A positive `fmax` switches each PSD to
625 : * its asymptotic power-law tail above that frequency; zero selects the multilayer default cutoff. Calculation and
626 : * write failures are collected by spatial-mode index and the first failure in grid order is returned after the
627 : * parallel loop. Files completed before a calculation or write failure are retained.
628 : *
629 : * \returns `error_t::noerror` when the complete grid is written, or a typed validation, calculation, or output
630 : * error.
631 : *
632 : */
633 : error_t makePSDGrid( const std::string &dir, ///< [in] the directory for output of the PSDs
634 : int mnMax, ///< [in] the positive maximum value of m and n in the grid
635 : realT dFreq, ///< [in] the positive temporal frequency spacing
636 : realT maxFreq, ///< [in] the positive maximum temporal frequency to calculate
637 : realT fmax = 0 ///< [in] maximum exactly calculated frequency, or 0 for the default cutoff
638 : );
639 :
640 : /// Analyze a PSD grid under closed-loop control.
641 : /** This always analyzes the simple integrator, and can also analyze the linear predictor controller. Outputs maps
642 : * of optimum gains, predictor coefficients, variances, and contrasts for a range of guide star magnitudes.
643 : * Optionally calculates speckle lifetimes. Optionally writes the closed-loop PSDs and transfer functions.
644 : */
645 : int analyzePSDGrid( const std::string &subDir, /**< [out] the sub-directory of psdDir where to write the
646 : results. Is created. */
647 : const std::string &psdDir, ///< [in] the directory containing the grid of PSDs.
648 : int mnMax, ///< [in] the maximum value of m and n in the grid.
649 : int mnCon, ///< [in] the maximum value of m and n which can be controlled.
650 : realT gfixed, ///< [in] if \> 0 then this fixed gain is used in the SI.
651 : int lpNc, /**< [in] the number of linear predictor coefficients to analyze.
652 : If 0 then LP is not analyzed.*/
653 : realT lpRegPrecision, /**< [in] the initial precision for the LP regularization
654 : algorithm. Normal value is 2. Higher is faster. Decrease
655 : if getting stuck in local minima.*/
656 : std::vector<realT> &mags, ///< [in] the guide star magnitudes to analyze for.
657 : int lifetimeTrials = 0, /**< [in] [optional] number of trials used for calculating speckle
658 : lifetimes. If 0,lifetimes are not calculated. */
659 : bool ucLifeTs = false, /**< [in] [optional] flag controlling whether lifetimes are
660 : calculated for uncontrolled modes.*/
661 : bool writePSDs = false, ///< [in] [optional] flag controlling if resultant PSDs are saved
662 : bool writeXfer = false /**< [in] [optional] flag controlling if resultant
663 : transfer functions are saved*/
664 : );
665 :
666 : int intensityPSD( const std::string &subDir, // sub-directory of psdDir which contains the controlled system
667 : // results, and where the lifetimes will be written.
668 : const std::string &psdDir, // directory containing the grid of PSDs
669 : const std::string &CvdPath, // path to the covariance decomposition
670 : int mnMax, ///< [in] the maximum value of m and n in the grid.
671 : int mnCon, ///< [in] the maximum value of m and n which can be controlled.
672 : std::vector<realT> &mags, ///< [in] the guide star magnitudes
673 : int lifetimeTrials, /**< [in] [optional] number of trials used for calculating
674 : speckle lifetimes. If 0, lifetimes are not
675 : calculated.*/
676 : bool writePSDs /**< [in] [optional] flag controlling if resultant
677 : PSDs are saved*/ );
678 :
679 : /** \name Disk Storage
680 : * These methods handle writing to and reading from disk. The calculated PSDs are store in the mx::BinVector binary
681 : * format.
682 : *
683 : * A grid of PSDs is specified by its directory name. The directory contains one frequency file (freq.binv), and a
684 : * set of PSD files, named according to psd_\<m\>_\<n\>_.binv.
685 : *
686 : *
687 : * @{
688 : */
689 : /// Get the frequency scale for a PSD grid.
690 : /**
691 : */
692 : int getGridFreq( std::vector<realT> &freq, ///< [out] the vector to populate with the frequency scale.
693 : const std::string &dir ///< [in] specifies the directory containing the grid.
694 : );
695 :
696 : /// Get a single PSD from a PSD grid.
697 : /**
698 : */
699 : int getGridPSD( std::vector<realT> &psd, ///< [out] the vector to populate with the PSD.
700 : const std::string &dir, ///< [in] specifies the directory containing the grid.
701 : int m, ///< [in] specifies the u component of spatial frequency.
702 : int n ///< [in] specifies the v component of spatial frequency.
703 : );
704 :
705 : /// Get both the frequency scale and a single PSD from a PSD grid.
706 : /**
707 : */
708 : int getGridPSD( std::vector<realT> &freq, ///< [out] the vector to populate with the frequency scale.
709 : std::vector<realT> &psd, ///< [out] the vector to populate with the PSD.
710 : const std::string &dir, ///< [in] specifies the directory containing the grid.
711 : int m, ///< [in] specifies the u component of spatial frequency.
712 : int n ///< [in] specifies the v component of spatial frequency.
713 : );
714 :
715 : ///@}
716 : };
717 :
718 : template <typename realT, typename aosysT>
719 34 : fourierTemporalPSD<realT, aosysT>::fourierTemporalPSD()
720 : {
721 34 : m_aosys = nullptr;
722 34 : initialize();
723 34 : }
724 :
725 : template <typename realT, typename aosysT>
726 13 : fourierTemporalPSD<realT, aosysT>::fourierTemporalPSD( fourierTemporalPSD_detail::gslWorkspaceAllocator allocator )
727 13 : : m_workspaceAllocator( allocator )
728 : {
729 13 : m_aosys = nullptr;
730 13 : initialize();
731 13 : }
732 :
733 : template <typename realT, typename aosysT>
734 12 : error_t fourierTemporalPSD<realT, aosysT>::allocateWorkspace()
735 : {
736 12 : if( m_workspace != nullptr )
737 : {
738 0 : return error_t::noerror;
739 : }
740 :
741 12 : if( m_workspaceAllocator == nullptr )
742 : {
743 0 : return internal::mxlib_error_report( error_t::invalidconfig, "GSL workspace allocator is null" );
744 : }
745 :
746 12 : m_workspace.reset( m_workspaceAllocator( WSZ ) );
747 12 : if( m_workspace == nullptr )
748 : {
749 15 : return internal::mxlib_error_report( error_t::allocerr, "could not allocate GSL integration workspace" );
750 : }
751 :
752 7 : return error_t::noerror;
753 : }
754 :
755 : template <typename realT, typename aosysT>
756 47 : void fourierTemporalPSD<realT, aosysT>::initialize()
757 : {
758 47 : _useBasis = basis::modified;
759 :
760 47 : _absTol = 1e-10;
761 47 : _relTol = 1e-4;
762 47 : }
763 :
764 : template <typename realT, typename aosysT>
765 21 : error_t fourierTemporalPSD<realT, aosysT>::validateAtmosphere( int layer_i )
766 : {
767 21 : if( m_aosys == nullptr )
768 : {
769 6 : return internal::mxlib_error_report( error_t::invalidconfig, "AO system pointer is null" );
770 : }
771 :
772 19 : if( !math::isFinite( m_aosys->D() ) || m_aosys->D() <= 0 )
773 : {
774 2 : return internal::mxlib_error_report( error_t::invalidconfig,
775 3 : "AO system aperture diameter must be finite and positive" );
776 : }
777 :
778 18 : auto &atmosphere = m_aosys->atm;
779 18 : const error_t atmosphereStatus = atmosphere.validate();
780 18 : if( atmosphereStatus != error_t::noerror )
781 : {
782 2 : return atmosphereStatus;
783 : }
784 :
785 16 : const size_t layerCount = atmosphere.n_layers();
786 42 : for( size_t index = 0; index < layerCount; ++index )
787 : {
788 27 : if( atmosphere.layer_v_wind( static_cast<int>( index ) ) <= 0 )
789 : {
790 2 : return internal::mxlib_error_report( error_t::invalidconfig,
791 3 : "Fourier temporal PSD layers require positive wind speed" );
792 : }
793 : }
794 :
795 15 : if( layer_i < -1 || ( layer_i >= 0 && static_cast<size_t>( layer_i ) >= layerCount ) )
796 : {
797 3 : return internal::mxlib_error_report( error_t::invalidarg, "atmosphere layer index is out of range" );
798 : }
799 :
800 14 : return error_t::noerror;
801 : }
802 :
803 : template <typename realT, typename aosysT>
804 32 : error_t fourierTemporalPSD<realT, aosysT>::validatePsdInputs( const std::vector<realT> &PSD,
805 : const std::vector<realT> &freq,
806 : realT m,
807 : realT n,
808 : int p,
809 : realT fmax,
810 : int layer_i,
811 : fourierTemporalPSDPolicy policy )
812 : {
813 32 : if( freq.empty() || PSD.size() != freq.size() )
814 : {
815 4 : return internal::mxlib_error_report( error_t::sizeerr,
816 6 : "PSD and frequency vectors must have the same nonzero size" );
817 : }
818 :
819 285 : for( size_t index = 0; index < freq.size(); ++index )
820 : {
821 258 : if( !math::isFinite( freq[index] ) || freq[index] < 0 || ( index > 0 && freq[index] <= freq[index - 1] ) )
822 : {
823 3 : return internal::mxlib_error_report(
824 3 : error_t::invalidarg,
825 9 : "frequency grid must be finite, nonnegative, and strictly increasing" );
826 : }
827 : }
828 :
829 27 : if( !math::isFinite( m ) || !math::isFinite( n ) || !math::isFinite( fmax ) || fmax < 0 )
830 : {
831 2 : return internal::mxlib_error_report(
832 2 : error_t::invalidarg,
833 6 : "mode coordinates must be finite and frequency cutoff must be finite and nonnegative" );
834 : }
835 :
836 25 : if( p != -1 && p != 1 )
837 : {
838 3 : return internal::mxlib_error_report( error_t::invalidarg, "Fourier-mode parity must be -1 or +1" );
839 : }
840 :
841 24 : if( _useBasis != basis::basic && _useBasis != basis::modified )
842 : {
843 3 : return internal::mxlib_error_report( error_t::invalidarg, "value of _useBasis is not valid" );
844 : }
845 :
846 23 : if( policy != fourierTemporalPSDPolicy::permissive && policy != fourierTemporalPSDPolicy::strict )
847 : {
848 0 : return internal::mxlib_error_report( error_t::invalidarg, "quadrature policy is not valid" );
849 : }
850 :
851 23 : if( !math::isFinite( _absTol ) || _absTol <= 0 || !math::isFinite( _relTol ) || _relTol <= 0 || _relTol >= 1 )
852 : {
853 2 : return internal::mxlib_error_report(
854 2 : error_t::invalidconfig,
855 6 : "GSL absolute tolerance must be positive and relative tolerance must be between zero and one" );
856 : }
857 :
858 21 : if( !math::isFinite( m_f0 ) || m_f0 < 0 )
859 : {
860 0 : return internal::mxlib_error_report( error_t::invalidconfig,
861 0 : "turbulence boiling parameter must be finite and nonnegative" );
862 : }
863 :
864 21 : const error_t atmosphereStatus = validateAtmosphere( layer_i );
865 21 : if( atmosphereStatus != error_t::noerror )
866 : {
867 7 : return atmosphereStatus;
868 : }
869 :
870 14 : if( _useBasis == basis::modified )
871 : {
872 0 : if( !math::isFinite( m_aosys->lam_sci() ) || m_aosys->lam_sci() <= 0 || !math::isFinite( m_aosys->lam_wfs() ) ||
873 0 : m_aosys->lam_wfs() <= 0 || !math::isFinite( m_aosys->zeta() ) ||
874 0 : std::abs( m_aosys->zeta() ) >= math::half_pi<realT>() )
875 : {
876 0 : return internal::mxlib_error_report(
877 0 : error_t::invalidconfig,
878 0 : "AO wavelengths must be positive and zenith angle must lie strictly between -pi/2 and pi/2" );
879 : }
880 : }
881 :
882 28 : if( !math::isFinite( m_aosys->spatialFilter_ku() ) || m_aosys->spatialFilter_ku() <= 0 ||
883 28 : !math::isFinite( m_aosys->spatialFilter_kv() ) || m_aosys->spatialFilter_kv() <= 0 )
884 : {
885 0 : return internal::mxlib_error_report( error_t::invalidconfig,
886 0 : "AO spatial-filter limits must be finite and positive" );
887 : }
888 :
889 14 : return error_t::noerror;
890 : }
891 :
892 : template <typename realT, typename aosysT>
893 1 : void fourierTemporalPSD<realT, aosysT>::absTol( realT at )
894 : {
895 1 : _absTol = at;
896 1 : }
897 :
898 : template <typename realT, typename aosysT>
899 0 : realT fourierTemporalPSD<realT, aosysT>::absTol()
900 : {
901 0 : return _absTol;
902 : }
903 :
904 : template <typename realT, typename aosysT>
905 1 : void fourierTemporalPSD<realT, aosysT>::relTol( realT rt )
906 : {
907 1 : _relTol = rt;
908 1 : }
909 :
910 : template <typename realT, typename aosysT>
911 0 : realT fourierTemporalPSD<realT, aosysT>::relTol()
912 : {
913 0 : return _relTol;
914 : }
915 :
916 : template <typename realT, typename aosysT>
917 4 : realT fourierTemporalPSD<realT, aosysT>::fastestPeak( int m, int n )
918 : {
919 : realT ku, kv, vu, vv;
920 :
921 4 : ku = ( (realT)m / m_aosys->D() );
922 4 : kv = ( (realT)n / m_aosys->D() );
923 :
924 4 : realT f, fmax = 0;
925 :
926 8 : for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
927 : {
928 4 : vu = m_aosys->atm.layer_v_wind( i ) * cos( m_aosys->atm.layer_dir( i ) );
929 4 : vv = m_aosys->atm.layer_v_wind( i ) * sin( m_aosys->atm.layer_dir( i ) );
930 :
931 4 : f = fabs( ku * vu + kv * vv );
932 4 : if( f > fmax )
933 3 : fmax = f;
934 : }
935 :
936 4 : return fmax;
937 : }
938 :
939 : template <typename realT, typename aosysT>
940 22 : error_t fourierTemporalPSD<realT, aosysT>::singleLayerPSD( std::vector<realT> &PSD,
941 : std::vector<realT> &freq,
942 : realT m,
943 : realT n,
944 : int layer_i,
945 : int p,
946 : realT fmax,
947 : reportT *report,
948 : fourierTemporalPSDPolicy policy )
949 : {
950 22 : reportT localReport;
951 22 : reportT &activeReport = report == nullptr ? localReport : *report;
952 22 : activeReport.clear();
953 :
954 22 : const error_t status = validatePsdInputs( PSD, freq, m, n, p, fmax, layer_i, policy );
955 22 : if( status != error_t::noerror )
956 : {
957 14 : return status;
958 : }
959 :
960 8 : fourierTemporalPSD_detail::scopedGslErrorHandlerOff handlerGuard;
961 8 : return singleLayerPSDImpl( PSD, freq, m, n, layer_i, p, fmax, activeReport, policy );
962 22 : }
963 :
964 : template <typename realT, typename aosysT>
965 12 : error_t fourierTemporalPSD<realT, aosysT>::singleLayerPSDImpl( std::vector<realT> &PSD,
966 : std::vector<realT> &freq,
967 : realT m,
968 : realT n,
969 : int layer_i,
970 : int p,
971 : realT fmax,
972 : reportT &report,
973 : fourierTemporalPSDPolicy policy )
974 : {
975 12 : if( fmax == 0 )
976 1 : fmax = freq[freq.size() - 1];
977 :
978 12 : if( freq[0] > fmax )
979 : {
980 1 : return internal::mxlib_error_report(
981 1 : error_t::invalidarg,
982 3 : "at least one exact frequency bin is required to initialize the PSD tail" );
983 : }
984 :
985 11 : realT v_wind = m_aosys->atm.layer_v_wind( layer_i );
986 11 : realT q_wind = m_aosys->atm.layer_dir( layer_i );
987 :
988 : // Rotate the basis
989 11 : realT cq = cos( q_wind );
990 11 : realT sq = sin( q_wind );
991 :
992 11 : realT scale = 2 * ( 1 / v_wind ); // Factor of 2 for negative frequencies.
993 :
994 : // Create a local instance so that we're reentrant
995 11 : fourierTemporalPSD<realT, aosysT> params( m_workspaceAllocator );
996 11 : const error_t workspaceStatus = params.allocateWorkspace();
997 11 : if( workspaceStatus != error_t::noerror )
998 : {
999 5 : return workspaceStatus;
1000 : }
1001 :
1002 6 : params.m_aosys = m_aosys;
1003 6 : params._layer_i = layer_i;
1004 6 : params.m_m = m * cq + n * sq;
1005 6 : params.m_n = -m * sq + n * cq;
1006 6 : params.m_cq = cq; // for de-rotating ku and kv for spatial filtering
1007 6 : params.m_sq = sq; // for de-rotation ku and kv for spatial filtering
1008 12 : if( m_aosys->spatialFilter_ku() < std::numeric_limits<realT>::max() ||
1009 6 : m_aosys->spatialFilter_kv() < std::numeric_limits<realT>::max() )
1010 0 : params.m_spatialFilter = true;
1011 :
1012 6 : params.m_p = p;
1013 6 : params.m_f0 = m_f0;
1014 :
1015 6 : params.m_mode_i = m_mode_i;
1016 6 : params.m_modeCoeffs = m_modeCoeffs;
1017 6 : params.m_minCoeffVal = m_minCoeffVal;
1018 :
1019 6 : realT result{ 0 };
1020 6 : realT error{ 0 };
1021 6 : error_t returnStatus = error_t::noerror;
1022 :
1023 : // Setup the GSL calculation
1024 : gsl_function func;
1025 6 : switch( _useBasis )
1026 : {
1027 6 : case basis::basic: // MXAO_FTPSD_BASIS_BASIC:
1028 6 : func.function = &F_basic<realT, aosysT>;
1029 6 : break;
1030 0 : case basis::modified: // MXAO_FTPSD_BASIS_MODIFIED:
1031 0 : func.function = &F_mod<realT, aosysT>;
1032 0 : break;
1033 0 : default:
1034 0 : return internal::mxlib_error_report( error_t::invalidarg, "value of _useBasis is not valid." );
1035 : }
1036 :
1037 6 : func.params = ¶ms;
1038 :
1039 : // Here we only calculate up to fmax.
1040 6 : size_t i = 0;
1041 106 : while( freq[i] <= fmax )
1042 : {
1043 103 : params.m_f = freq[i];
1044 :
1045 103 : const int ec = gsl_integration_qagi( &func, _absTol, _relTol, WSZ, params.m_workspace.get(), &result, &error );
1046 103 : report.record( ec, static_cast<size_t>( layer_i ), freq[i], result, error, _absTol, _relTol );
1047 :
1048 103 : const error_t integrationStatus = fourierTemporalPSD_detail::applyPolicy( ec, policy );
1049 103 : if( integrationStatus != error_t::noerror )
1050 : {
1051 0 : if( !fourierTemporalPSD_detail::isConvergenceStatus( ec ) )
1052 : {
1053 0 : return internal::mxlib_error_report( integrationStatus,
1054 0 : std::string( "gsl_integration_qagi failed: " ) +
1055 0 : gsl_strerror( ec ) );
1056 : }
1057 :
1058 0 : returnStatus = integrationStatus;
1059 : }
1060 :
1061 103 : if( !math::isFinite( result ) || !math::isFinite( error ) )
1062 : {
1063 0 : return internal::mxlib_error_report( error_t::liberr,
1064 0 : "gsl_integration_qagi returned a nonfinite result or error estimate" );
1065 : }
1066 :
1067 103 : PSD[i] = scale * result;
1068 :
1069 103 : ++i;
1070 103 : if( i >= freq.size() )
1071 3 : break;
1072 : }
1073 :
1074 : // Now fill in from fmax to the actual max frequency with a -(alpha+2) power law.
1075 6 : size_t j = i;
1076 :
1077 6 : if( j == freq.size() )
1078 3 : return returnStatus;
1079 :
1080 : // First average up to the last 50 exactly integrated bins after projecting them to the first tail frequency.
1081 3 : constexpr size_t maximumTailAverageCount = 50;
1082 3 : const size_t tailAverageCount = std::min( i, maximumTailAverageCount );
1083 3 : PSD[j] = 0;
1084 103 : for( size_t k = tailAverageCount; k > 0; --k )
1085 : {
1086 200 : PSD[j] +=
1087 100 : PSD[i - k] * pow( freq[i - k] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds<realT>());
1088 : }
1089 3 : PSD[j] /= static_cast<realT>( tailAverageCount );
1090 3 : ++j;
1091 3 : ++i;
1092 3 : if( j == freq.size() )
1093 0 : return returnStatus;
1094 65 : while( j < freq.size() )
1095 : {
1096 124 : PSD[j] =
1097 62 : PSD[i - 1] * pow( freq[i - 1] / freq[j], m_aosys->atm.alpha( layer_i ) + 2 ); // seventeen_thirds<realT>());
1098 62 : ++j;
1099 : }
1100 :
1101 3 : return returnStatus;
1102 11 : }
1103 :
1104 : template <typename realT, typename aosysT>
1105 0 : error_t fourierTemporalPSD<realT, aosysT>::multiLayerPSDImpl( std::vector<realT> &PSD,
1106 : std::vector<realT> &freq,
1107 : realT m,
1108 : realT n,
1109 : int p,
1110 : realT fmax,
1111 : reportT &report,
1112 : fourierTemporalPSDPolicy policy,
1113 : isParallel<true> parallel )
1114 : {
1115 : static_cast<void>( parallel );
1116 :
1117 0 : const size_t layerCount = m_aosys->atm.n_layers();
1118 0 : std::vector<error_t> layerStatus( layerCount, error_t::noerror );
1119 0 : std::vector<reportT> layerReport( layerCount );
1120 :
1121 0 : #pragma omp parallel
1122 : {
1123 : // Records each layer PSD
1124 : std::vector<realT> single_PSD( freq.size() );
1125 :
1126 : #pragma omp for
1127 : for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
1128 : {
1129 : std::fill( single_PSD.begin(), single_PSD.end(), 0 );
1130 : layerStatus[i] =
1131 : singleLayerPSDImpl( single_PSD, freq, m, n, static_cast<int>( i ), p, fmax, layerReport[i], policy );
1132 :
1133 : // Now add the single layer PSD to the overall PSD, weighted by Cn2
1134 : #pragma omp critical
1135 : if( layerStatus[i] == error_t::noerror )
1136 : {
1137 : for( size_t j = 0; j < freq.size(); ++j )
1138 : {
1139 : PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j];
1140 : }
1141 : }
1142 : }
1143 : }
1144 :
1145 0 : error_t returnStatus = error_t::noerror;
1146 0 : for( size_t i = 0; i < layerCount; ++i )
1147 : {
1148 0 : report.merge( layerReport[i] );
1149 0 : if( returnStatus == error_t::noerror && layerStatus[i] != error_t::noerror )
1150 : {
1151 0 : returnStatus = layerStatus[i];
1152 : }
1153 : }
1154 :
1155 0 : return returnStatus;
1156 0 : }
1157 :
1158 : template <typename realT, typename aosysT>
1159 4 : error_t fourierTemporalPSD<realT, aosysT>::multiLayerPSDImpl( std::vector<realT> &PSD,
1160 : std::vector<realT> &freq,
1161 : realT m,
1162 : realT n,
1163 : int p,
1164 : realT fmax,
1165 : reportT &report,
1166 : fourierTemporalPSDPolicy policy,
1167 : isParallel<false> parallel )
1168 : {
1169 : static_cast<void>( parallel );
1170 :
1171 : // Records each layer PSD
1172 4 : std::vector<realT> single_PSD( freq.size() );
1173 4 : error_t returnStatus = error_t::noerror;
1174 :
1175 8 : for( size_t i = 0; i < m_aosys->atm.n_layers(); ++i )
1176 : {
1177 4 : std::fill( single_PSD.begin(), single_PSD.end(), 0 );
1178 4 : reportT layerReport;
1179 : const error_t layerStatus =
1180 4 : singleLayerPSDImpl( single_PSD, freq, m, n, static_cast<int>( i ), p, fmax, layerReport, policy );
1181 4 : report.merge( layerReport );
1182 4 : if( returnStatus == error_t::noerror && layerStatus != error_t::noerror )
1183 : {
1184 4 : returnStatus = layerStatus;
1185 : }
1186 :
1187 : // Now add the single layer PSD to the overall PSD, weighted by Cn2
1188 4 : if( layerStatus == error_t::noerror )
1189 : {
1190 0 : for( size_t j = 0; j < freq.size(); ++j )
1191 : {
1192 0 : PSD[j] += m_aosys->atm.layer_Cn2( i ) * single_PSD[j];
1193 : }
1194 : }
1195 : }
1196 :
1197 4 : return returnStatus;
1198 4 : }
1199 :
1200 : template <typename realT, typename aosysT>
1201 : template <bool parallel>
1202 7 : error_t fourierTemporalPSD<realT, aosysT>::multiLayerPSD( std::vector<realT> &PSD,
1203 : std::vector<realT> &freq,
1204 : realT m,
1205 : realT n,
1206 : int p,
1207 : realT fmax,
1208 : reportT *report,
1209 : fourierTemporalPSDPolicy policy )
1210 : {
1211 7 : reportT localReport;
1212 7 : reportT &activeReport = report == nullptr ? localReport : *report;
1213 7 : activeReport.clear();
1214 :
1215 7 : const error_t validationStatus = validatePsdInputs( PSD, freq, m, n, p, fmax, -1, policy );
1216 7 : if( validationStatus != error_t::noerror )
1217 : {
1218 3 : return validationStatus;
1219 : }
1220 :
1221 : // PSD is zeroed every time to make sure we don't accumulate on repeated calls
1222 8 : for( size_t j = 0; j < PSD.size(); ++j )
1223 4 : PSD[j] = 0;
1224 :
1225 4 : if( fmax == 0 )
1226 : {
1227 4 : fmax = 150 + 2 * fastestPeak( m, n );
1228 : }
1229 :
1230 4 : fourierTemporalPSD_detail::scopedGslErrorHandlerOff handlerGuard;
1231 4 : return multiLayerPSDImpl( PSD, freq, m, n, p, fmax, activeReport, policy, isParallel<parallel>() );
1232 7 : }
1233 :
1234 : template <typename realT, typename aosysT>
1235 10 : error_t fourierTemporalPSD<realT, aosysT>::makePSDGrid(
1236 : const std::string &dir, int mnMax, realT dFreq, realT maxFreq, realT fmax )
1237 : {
1238 10 : if( dir.empty() )
1239 : {
1240 3 : return internal::mxlib_error_report( error_t::invalidarg, "PSD grid output directory must not be empty" );
1241 : }
1242 :
1243 8 : if( mnMax <= 0 || !math::isFinite( dFreq ) || dFreq <= 0 || !math::isFinite( maxFreq ) || maxFreq <= 0 ||
1244 17 : !math::isFinite( fmax ) || fmax < 0 )
1245 : {
1246 5 : return internal::mxlib_error_report(
1247 5 : error_t::invalidarg,
1248 15 : "PSD grid extent and frequency controls must be finite and positive, with a nonnegative cutoff" );
1249 : }
1250 :
1251 4 : const realT sampleCount = maxFreq / dFreq;
1252 4 : if( !math::isFinite( sampleCount ) || sampleCount > static_cast<realT>( std::numeric_limits<int>::max() ) )
1253 : {
1254 3 : return internal::mxlib_error_report( error_t::sizeerr, "PSD grid sample count exceeds the supported range" );
1255 : }
1256 :
1257 6 : const std::vector<realT> validationFrequency{ 0 };
1258 6 : const std::vector<realT> validationPsd{ 0 };
1259 3 : const error_t validationStatus = validatePsdInputs( validationPsd,
1260 : validationFrequency,
1261 : 0,
1262 : 0,
1263 : 1,
1264 : fmax,
1265 : -1,
1266 : fourierTemporalPSDPolicy::permissive );
1267 3 : if( validationStatus != error_t::noerror )
1268 : {
1269 1 : return validationStatus;
1270 : }
1271 :
1272 2 : std::vector<realT> freq;
1273 :
1274 2 : std::vector<sigproc::fourierModeDef> spf;
1275 :
1276 2 : std::string fn;
1277 :
1278 2 : sigproc::makeFourierModeFreqs_Rect( spf, 2 * mnMax );
1279 :
1280 : // Calculate number of samples, and make sure we get to at least maxFreq
1281 2 : int N = (int)( maxFreq / dFreq );
1282 2 : if( N * dFreq < maxFreq )
1283 0 : N += 1;
1284 :
1285 : /*** Dump Params to file ***/
1286 2 : error_t status = ioutils::createDirectories( dir );
1287 2 : if( status != error_t::noerror )
1288 : {
1289 3 : return internal::mxlib_error_report( status, "could not create PSD grid output directory" );
1290 : }
1291 :
1292 1 : std::ofstream fout;
1293 1 : fn = dir + '/' + "params.txt";
1294 1 : fout.open( fn );
1295 1 : if( !fout.is_open() )
1296 : {
1297 0 : return internal::mxlib_error_report( error_t::fileoerr, "could not open PSD grid parameter file" );
1298 : }
1299 :
1300 1 : fout << "#---------------------------\n";
1301 1 : m_aosys->dumpAOSystem( fout );
1302 1 : fout << "#---------------------------\n";
1303 1 : fout << "# PSD Grid Parameters\n";
1304 1 : fout << "# absTol " << _absTol << '\n';
1305 1 : fout << "# relTol " << _relTol << '\n';
1306 1 : fout << "# useBasis " << _useBasis << '\n';
1307 1 : fout << "# makePSDGrid call:\n";
1308 1 : fout << "# mnMax = " << mnMax << '\n';
1309 1 : fout << "# dFreq = " << dFreq << '\n';
1310 1 : fout << "# maxFreq = " << maxFreq << '\n';
1311 1 : fout << "# fmax = " << fmax << '\n';
1312 1 : fout << "#---------------------------\n";
1313 :
1314 1 : fout.close();
1315 1 : if( !fout )
1316 : {
1317 0 : return internal::mxlib_error_report( error_t::filewerr, "could not write PSD grid parameter file" );
1318 : }
1319 :
1320 : // Make directory
1321 1 : std::string psddir = dir + "/psds";
1322 1 : status = ioutils::createDirectories( psddir );
1323 1 : if( status != error_t::noerror )
1324 : {
1325 0 : return internal::mxlib_error_report( status, "could not create PSD output directory" );
1326 : }
1327 :
1328 : // Create frequency scale.
1329 1 : math::vectorScale( freq, N, dFreq, 0 ); // dFreq); //offset from 0 by dFreq, so f=0 not included
1330 :
1331 1 : fn = psddir + '/' + "freq.binv";
1332 :
1333 1 : if( ioutils::writeBinVector( fn, freq ) != 0 )
1334 : {
1335 0 : return internal::mxlib_error_report( error_t::filewerr, "could not write PSD frequency grid" );
1336 : }
1337 :
1338 1 : size_t nLoops = 0.5 * spf.size();
1339 1 : std::vector<error_t> modeStatus( nLoops, error_t::noerror );
1340 :
1341 1 : ipc::ompLoopWatcher<> watcher( nLoops, std::cout );
1342 :
1343 1 : #pragma omp parallel
1344 : {
1345 : std::vector<realT> PSD;
1346 : PSD.resize( N );
1347 : std::string fname;
1348 :
1349 : int m, n;
1350 :
1351 : #pragma omp for
1352 : for( size_t i = 0; i < nLoops; ++i )
1353 : {
1354 : m = spf[i * 2].m;
1355 : n = spf[i * 2].n;
1356 :
1357 : if( fabs( (realT)m / m_aosys->D() ) >= m_aosys->spatialFilter_ku() ||
1358 : fabs( (realT)n / m_aosys->D() ) >= m_aosys->spatialFilter_kv() )
1359 : {
1360 : watcher.incrementAndOutputStatus();
1361 : continue;
1362 : }
1363 :
1364 : modeStatus[i] = multiLayerPSD<false>( PSD, freq, m, n, 1, fmax );
1365 : if( modeStatus[i] != error_t::noerror )
1366 : {
1367 : watcher.incrementAndOutputStatus();
1368 : continue;
1369 : }
1370 :
1371 : fname = std::format( "{}/psd_{}_{}.binv", psddir, m, n );
1372 : // psddir + '/' + "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) +
1373 : // ".binv";
1374 :
1375 : if( ioutils::writeBinVector( fname, PSD ) != 0 )
1376 : {
1377 : modeStatus[i] = error_t::filewerr;
1378 : }
1379 :
1380 : watcher.incrementAndOutputStatus();
1381 : }
1382 : }
1383 :
1384 1 : for( size_t index = 0; index < modeStatus.size(); ++index )
1385 : {
1386 1 : if( modeStatus[index] != error_t::noerror )
1387 : {
1388 1 : return modeStatus[index];
1389 : }
1390 : }
1391 :
1392 0 : return error_t::noerror;
1393 3 : }
1394 :
1395 : template <typename realT, typename aosysT>
1396 0 : int fourierTemporalPSD<realT, aosysT>::analyzePSDGrid( const std::string &subDir,
1397 : const std::string &psdDir,
1398 : int mnMax,
1399 : int mnCon,
1400 : realT gfixed,
1401 : int lpNc,
1402 : realT lpRegPrecision,
1403 : std::vector<realT> &mags,
1404 : int lifetimeTrials,
1405 : bool uncontrolledLifetimes,
1406 : bool writePSDs,
1407 : bool writeXfer )
1408 : {
1409 :
1410 0 : std::string dir = psdDir + "/" + subDir;
1411 :
1412 : /*** Dump Params to file ***/
1413 0 : mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
1414 :
1415 0 : std::ofstream fout;
1416 0 : std::string fn = dir + '/' + "params.txt";
1417 0 : fout.open( fn );
1418 :
1419 0 : fout << "#---------------------------\n";
1420 0 : m_aosys->dumpAOSystem( fout );
1421 0 : fout << "#---------------------------\n";
1422 0 : fout << "# Analysis Parameters\n";
1423 0 : fout << "# mnMax = " << mnMax << '\n';
1424 0 : fout << "# mnCon = " << mnCon << '\n';
1425 0 : fout << "# lpNc = " << lpNc << '\n';
1426 0 : fout << "# mags = ";
1427 0 : for( size_t i = 0; i < mags.size() - 1; ++i )
1428 0 : fout << mags[i] << ",";
1429 0 : fout << mags[mags.size() - 1] << '\n';
1430 0 : fout << "# lifetimeTrials = " << lifetimeTrials << '\n';
1431 0 : fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n';
1432 0 : fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n';
1433 0 : fout << "# writeXfer = " << std::boolalpha << writeXfer << '\n';
1434 :
1435 0 : fout.close();
1436 :
1437 : //**** Calculating A Variance Map ****//
1438 :
1439 0 : realT fs = 1.0 / m_aosys->tauWFS();
1440 0 : realT tauWFS = m_aosys->tauWFS();
1441 0 : realT deltaTau = m_aosys->deltaTau();
1442 :
1443 0 : std::vector<sigproc::fourierModeDef> fms;
1444 :
1445 0 : sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax );
1446 0 : size_t nModes = 0.5 * fms.size();
1447 :
1448 0 : Eigen::Array<realT, -1, -1> gains, vars, speckleLifetimes, gains_lp, vars_lp, speckleLifetimes_lp;
1449 :
1450 0 : gains.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1451 0 : vars.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1452 0 : speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1453 :
1454 0 : gains( mnMax, mnMax ) = 0;
1455 0 : vars( mnMax, mnMax ) = 0;
1456 0 : speckleLifetimes( mnMax, mnMax ) = 0;
1457 :
1458 0 : gains_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1459 0 : vars_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1460 0 : speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
1461 :
1462 0 : gains_lp( mnMax, mnMax ) = 0;
1463 0 : vars_lp( mnMax, mnMax ) = 0;
1464 0 : speckleLifetimes_lp( mnMax, mnMax ) = 0;
1465 :
1466 0 : bool doLP = false;
1467 0 : if( lpNc > 1 )
1468 0 : doLP = true;
1469 0 : Eigen::Array<realT, -1, -1> lpC;
1470 :
1471 0 : if( doLP )
1472 : {
1473 0 : lpC.resize( nModes, lpNc );
1474 0 : lpC.setZero();
1475 : }
1476 :
1477 0 : std::vector<realT> S_si, S_lp;
1478 :
1479 0 : if( writePSDs )
1480 : {
1481 0 : for( size_t s = 0; s < mags.size(); ++s )
1482 : {
1483 0 : std::string psdOutDir = std::format( "{}/outputPSDS_{}_si", dir, mags[s] );
1484 : // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si";
1485 0 : mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
1486 :
1487 0 : if( doLP )
1488 : {
1489 0 : std::string psdOutDir = std::format( "{}/outputPSDS_{}_lp", dir, mags[s] );
1490 : // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp";
1491 0 : mkdir( psdOutDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
1492 0 : }
1493 : }
1494 : }
1495 :
1496 0 : m_aosys->beta_p( 1, 1 );
1497 :
1498 0 : ipc::ompLoopWatcher<> watcher( nModes * mags.size(), std::cout );
1499 0 : std::atomic<int> analysisStatus{ static_cast<int>( error_t::noerror ) };
1500 :
1501 0 : for( size_t s = 0; s < mags.size(); ++s )
1502 : {
1503 0 : m_aosys->starMag( mags[s] );
1504 :
1505 : // In non-parallel space calculate OG=Strehl
1506 :
1507 0 : realT opticalGain{ 1.0 };
1508 :
1509 0 : if( m_strehlOG )
1510 : {
1511 : // Iterative optical gain
1512 : /// \todo need upstream NCP and CP NCP and NCP NCP
1513 0 : realT ncp = m_aosys->ncp_wfe(); // save ncp and then set it to zero for this part.
1514 0 : m_aosys->ncp_wfe( 0 );
1515 :
1516 0 : realT lam_sci = m_aosys->lam_sci();
1517 0 : m_aosys->lam_sci( m_aosys->lam_wfs() );
1518 :
1519 0 : realT S = m_aosys->strehl();
1520 0 : std::cerr << S << "\n";
1521 :
1522 0 : for( int s = 0; s < 4; ++s )
1523 : {
1524 0 : m_aosys->opticalGain( S );
1525 0 : m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
1526 0 : S = m_aosys->strehl();
1527 0 : std::cerr << S << "\n";
1528 : }
1529 :
1530 0 : opticalGain = S;
1531 :
1532 0 : m_aosys->lam_sci( lam_sci );
1533 0 : m_aosys->ncp_wfe( ncp );
1534 : }
1535 :
1536 0 : m_aosys->optd( m_aosys->optd() ); // just trigger a re-calc
1537 0 : realT strehl = m_aosys->strehl();
1538 :
1539 0 : #pragma omp parallel
1540 : {
1541 : realT localMag = mags[s];
1542 :
1543 : realT var0;
1544 :
1545 : realT gopt, var;
1546 :
1547 : realT gopt_lp, var_lp;
1548 :
1549 : std::vector<realT> tfreq; // The frequency scale of the PSDs
1550 : std::vector<realT> tPSDp; // The open-loop turbulence PSD for a Fourier mode
1551 : std::vector<realT>
1552 : tPSDpPOL; // The pseudo-open-loop turbulence PSD for a Fourier mode, with optical gain effects included
1553 :
1554 : std::vector<realT> tfreqHF; // The above-Nyquist frequencies, saved if outputing the PSDS.
1555 : std::vector<realT> tPSDpHF; // The above-Nyquist component of the open-loop PSD, saved if outputing the
1556 : // PSDs.
1557 :
1558 : //**< Get the frequency grid, and nyquist limit it to f_s/2
1559 : getGridPSD( tfreq, tPSDp, psdDir, 0, 1 ); // To get the freq grid
1560 :
1561 : size_t imax = 0;
1562 : while( tfreq[imax] <= 0.5 * fs )
1563 : {
1564 : ++imax;
1565 : if( imax > tfreq.size() - 1 )
1566 : break;
1567 : }
1568 :
1569 : if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
1570 : {
1571 : ++imax;
1572 : }
1573 :
1574 : if( writePSDs )
1575 : {
1576 : tfreqHF.assign( tfreq.begin(), tfreq.end() );
1577 : }
1578 :
1579 : tfreq.erase( tfreq.begin() + imax, tfreq.end() );
1580 :
1581 : tPSDpPOL.resize( tfreq.size() ); // pre=allocate
1582 : //**>
1583 :
1584 : std::vector<realT> tPSDn; // The open-loop WFS noise PSD
1585 : tPSDn.resize( tfreq.size() );
1586 :
1587 : //**< Setup the controllers
1588 : mx::AO::analysis::clAOLinearPredictor<realT> tflp;
1589 : tflp.m_precision0 = lpRegPrecision;
1590 : /*tflp.m_min_sc0 = 0;
1591 : tflp.m_max_sc0 = 1000;*/
1592 :
1593 : mx::AO::analysis::clGainOpt<realT> go_si( tauWFS, deltaTau );
1594 : mx::AO::analysis::clGainOpt<realT> go_lp( tauWFS, deltaTau );
1595 :
1596 : go_si.f( tfreq );
1597 : go_lp.f( tfreq );
1598 :
1599 : realT gmax = 0;
1600 : realT gmax_lp = 0;
1601 : //**>
1602 :
1603 : int m, n;
1604 :
1605 : //**< For use in lifetime calculations
1606 : sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
1607 : std::vector<std::complex<realT>> ETFxn;
1608 : std::vector<std::complex<realT>> NTFxn;
1609 :
1610 : if( lifetimeTrials > 0 )
1611 : {
1612 : ETFxn.resize( tfreq.size() );
1613 : NTFxn.resize( tfreq.size() );
1614 :
1615 : if( writeXfer )
1616 : {
1617 : std::string tfOutFile = std::format( "{}/outputTF_{}_si", dir, mags[s] );
1618 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
1619 : ioutils::createDirectories( tfOutFile );
1620 : }
1621 :
1622 : if( doLP )
1623 : {
1624 : if( writeXfer )
1625 : {
1626 : std::string tfOutFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] );
1627 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
1628 : ioutils::createDirectories( tfOutFile );
1629 : }
1630 : }
1631 : }
1632 :
1633 : //**>
1634 :
1635 : /*#pragma omp critical
1636 : std::cerr << __FILE__ << " " << __LINE__ << "\n";
1637 : */
1638 : // want to schedule dynamic with small chunks so maximal processor usage,
1639 : // otherwise we can end up with a small number of cores being used at the end
1640 : #pragma omp for schedule( dynamic, 5 )
1641 : for( size_t i = 0; i < nModes; ++i )
1642 : {
1643 : // Determine the spatial frequency at this step
1644 : m = fms[2 * i].m;
1645 : n = fms[2 * i].n;
1646 :
1647 : if( fabs( (realT)m / m_aosys->D() ) >= m_aosys->spatialFilter_ku() ||
1648 : fabs( (realT)n / m_aosys->D() ) >= m_aosys->spatialFilter_kv() )
1649 : {
1650 : gains( mnMax + m, mnMax + n ) = 0;
1651 : gains( mnMax - m, mnMax - n ) = 0;
1652 :
1653 : gains_lp( mnMax + m, mnMax + n ) = 0;
1654 : gains_lp( mnMax - m, mnMax - n ) = 0;
1655 :
1656 : vars( mnMax + m, mnMax + n ) = 0;
1657 : vars( mnMax - m, mnMax - n ) = 0;
1658 :
1659 : vars_lp( mnMax + m, mnMax + n ) = 0;
1660 : vars_lp( mnMax - m, mnMax - n ) = 0;
1661 : speckleLifetimes( mnMax + m, mnMax + n ) = 0;
1662 : speckleLifetimes( mnMax - m, mnMax - n ) = 0;
1663 : speckleLifetimes_lp( mnMax + m, mnMax + n ) = 0;
1664 : speckleLifetimes_lp( mnMax - m, mnMax - n ) = 0;
1665 : }
1666 : else
1667 : {
1668 :
1669 : realT k = sqrt( m * m + n * n ) / m_aosys->D();
1670 :
1671 : //**< Get the open-loop turb. PSD
1672 : getGridPSD( tPSDp, psdDir, m, n );
1673 :
1674 : // Get integral of entire open-loop PSD
1675 : var0 = sigproc::psdVar( tfreq, tPSDp );
1676 :
1677 : if( writePSDs )
1678 : {
1679 : tPSDpHF.assign( tPSDp.begin() + imax, tPSDp.end() );
1680 : }
1681 :
1682 : // erase points above Nyquist limit
1683 : tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() );
1684 :
1685 : // And now determine the variance which has been erased.
1686 : // limVar is the out-of-band variance, which we add back in for completeness
1687 : realT limVar = 0; // var0 - sigproc::psdVar( tfreq, tPSDp);
1688 :
1689 : // And construct the POL psd
1690 : if( m_uncorrectedOG )
1691 : {
1692 : for( size_t n = 0; n < tPSDp.size(); ++n )
1693 : {
1694 : tPSDpPOL[n] = tPSDp[n] * pow( opticalGain, 2 );
1695 : }
1696 : }
1697 : else
1698 : {
1699 : for( size_t n = 0; n < tPSDp.size(); ++n )
1700 : {
1701 : tPSDpPOL[n] = tPSDp[n];
1702 : }
1703 : }
1704 : //**>
1705 :
1706 : //**< Determine if we're inside the hardwarecontrol limit
1707 : bool inside = false;
1708 :
1709 : if( m_aosys->circularLimit() )
1710 : {
1711 : if( m * m + n * n <= mnCon * mnCon )
1712 : inside = true;
1713 : }
1714 : else
1715 : {
1716 : if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
1717 : inside = true;
1718 : }
1719 : //**>
1720 :
1721 : /* This is to select specific points for troubleshooting*/
1722 : // if( !( ( m ==-2 && n == 84 ) || (m==-2 && n == 2800)) ) inside = false;
1723 :
1724 : if( inside )
1725 : {
1726 : // Get the WFS noise PSD (which is already resized to match tfreq)
1727 : wfsNoisePSD<realT>( tPSDn,
1728 : m_aosys->beta_p( m, n ) / sqrt( opticalGain ),
1729 : m_aosys->Fg( localMag ),
1730 : tauWFS,
1731 : m_aosys->npix_wfs( (size_t)0 ),
1732 : m_aosys->Fbg( (size_t)0 ),
1733 : m_aosys->ron_wfs( (size_t)0 ) );
1734 :
1735 : gmax = 0;
1736 : if( gfixed > 0 )
1737 : {
1738 : gopt = gfixed;
1739 : var = go_si.clVariance( tPSDp, tPSDn, gopt );
1740 : }
1741 : else
1742 : {
1743 : // Calculate gain using the POL PSD
1744 : error_t gainStatus = go_si.optGainOpenLoop( gopt, var, tPSDpPOL, tPSDn, true );
1745 : if( gainStatus != error_t::noerror )
1746 : {
1747 : int expected = static_cast<int>( error_t::noerror );
1748 : analysisStatus.compare_exchange_strong( expected, static_cast<int>( gainStatus ) );
1749 : gopt = 0;
1750 : var = go_si.clVariance( tPSDp, tPSDn, gopt );
1751 : }
1752 :
1753 : if( m_uncorrectedOG )
1754 : {
1755 : gopt *= opticalGain;
1756 : }
1757 :
1758 : // But use the variance from the actual POL
1759 : var = go_si.clVariance( tPSDp, tPSDn, gopt );
1760 :
1761 : // Output gain curve for this mode (for troubleshooting)
1762 : /*#pragma omp critical
1763 : {
1764 : std::string foutn = "gcurve_";
1765 : foutn += std::to_string(m) + "_" + std::to_string(n) + ".dat";
1766 :
1767 : std::ofstream foutf(foutn);
1768 :
1769 : for(size_t n = 0; n < 10000; ++n)
1770 : {
1771 : realT gg = (1.0*n)/10000.;
1772 : foutf << gg << " " << go_si.clVariance(tPSDp, tPSDn, gg) << "\n";
1773 : }
1774 :
1775 : foutf.close();
1776 :
1777 : std::cerr << "\n" << gmax << " " << gopt << " " << var << " " << go_si.clVariance(tPSDp,
1778 : tPSDn, 0.64) << "\n";
1779 : }*/
1780 : }
1781 :
1782 : var += limVar;
1783 :
1784 : if( doLP )
1785 : {
1786 : realT min_sc;
1787 : error_t rv = tflp.regularizeCoefficients( gmax_lp,
1788 : gopt_lp,
1789 : var_lp,
1790 : min_sc,
1791 : go_lp,
1792 : tPSDpPOL,
1793 : tPSDn,
1794 : lpNc );
1795 :
1796 : if( rv != error_t::noerror )
1797 : {
1798 : int expected = static_cast<int>( error_t::noerror );
1799 : analysisStatus.compare_exchange_strong( expected, static_cast<int>( rv ) );
1800 : }
1801 :
1802 : for( int n = 0; n < lpNc; ++n )
1803 : {
1804 : lpC( i, n ) = go_lp.a( n );
1805 : }
1806 :
1807 : if( m_uncorrectedOG )
1808 : {
1809 : go_lp.aScale( opticalGain );
1810 : go_lp.bScale( opticalGain );
1811 : gopt_lp *= opticalGain;
1812 : }
1813 :
1814 : var_lp = go_lp.clVariance( tPSDp, tPSDn, gopt_lp );
1815 : var_lp += limVar;
1816 : }
1817 : else
1818 : {
1819 : gopt_lp = 0;
1820 : }
1821 : }
1822 : else
1823 : {
1824 : // Zero the noise PSD
1825 : tPSDn.assign( tPSDn.size(), 0.0 );
1826 :
1827 : gopt = 0;
1828 : var = var0;
1829 : var_lp = var0;
1830 : gopt_lp = 0;
1831 : go_lp.a( std::vector<realT>( { 1 } ) );
1832 : go_lp.b( std::vector<realT>( { 1 } ) );
1833 : }
1834 :
1835 : //**< Determine if closed-loop is making a difference:
1836 :
1837 : if( gopt > 0 && var > var0 )
1838 : {
1839 : gopt = 0;
1840 : var = var0;
1841 : }
1842 :
1843 : if( gopt_lp > gopt && var_lp > var )
1844 : {
1845 : // Set LP to SI (or off if SI is off)
1846 : gopt_lp = gopt;
1847 : var_lp = var;
1848 : go_lp.a( std::vector<realT>( { 1 } ) );
1849 : go_lp.b( std::vector<realT>( { 1 } ) );
1850 : }
1851 : //**>
1852 :
1853 : //**< Fill in the gain and variance maps
1854 : gains( mnMax + m, mnMax + n ) = gopt;
1855 : gains( mnMax - m, mnMax - n ) = gopt;
1856 :
1857 : gains_lp( mnMax + m, mnMax + n ) = gopt_lp;
1858 : gains_lp( mnMax - m, mnMax - n ) = gopt_lp;
1859 :
1860 : vars( mnMax + m, mnMax + n ) = var;
1861 : vars( mnMax - m, mnMax - n ) = var;
1862 :
1863 : vars_lp( mnMax + m, mnMax + n ) = var_lp;
1864 : vars_lp( mnMax - m, mnMax - n ) = var_lp;
1865 : //**>
1866 :
1867 : //**< Calculate Speckle Lifetimes
1868 : if( ( lifetimeTrials > 0 || writeXfer ) && ( uncontrolledLifetimes || inside ) )
1869 : {
1870 : std::vector<realT> spfreq, sppsd;
1871 :
1872 : if( gopt > 0 )
1873 : {
1874 : for( size_t i = 0; i < tfreq.size(); ++i )
1875 : {
1876 : ETFxn[i] = go_si.clETF( i, gopt );
1877 : NTFxn[i] = go_si.clNTF( i, gopt );
1878 : }
1879 : }
1880 : else
1881 : {
1882 : for( size_t i = 0; i < tfreq.size(); ++i )
1883 : {
1884 : ETFxn[i] = 1;
1885 : NTFxn[i] = 0;
1886 : }
1887 : }
1888 :
1889 : if( writeXfer )
1890 : {
1891 : std::string tfOutFile = std::format( "{}/outputTF_{}_si", dir, mags[s] );
1892 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
1893 :
1894 : std::string etfOutFile = std::format( "{}/etf_{}_{}.binv", tfOutFile, m, n );
1895 : // tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' +
1896 : // ioutils::convert ToString( n ) + ".binv";
1897 : ioutils::writeBinVector( etfOutFile, ETFxn );
1898 :
1899 : std::string ntfOutFile = std::format( "{}/ntf_{}_{}.binv", tfOutFile, m, n );
1900 : // tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' +
1901 : // ioutils::convert ToString( n ) + ".binv";
1902 : ioutils::writeBinVector( ntfOutFile, NTFxn );
1903 :
1904 : if( i == 0 ) // Write freq on the first one
1905 : {
1906 : std::string fOutFile = tfOutFile + "freq.binv";
1907 : ioutils::writeBinVector( fOutFile, tfreq );
1908 : }
1909 : }
1910 :
1911 : if( lifetimeTrials > 0 )
1912 : {
1913 : speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1914 : realT spvar = sigproc::psdVar( spfreq, sppsd );
1915 :
1916 : realT splifeT = 100.0;
1917 : realT error;
1918 :
1919 : realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1920 :
1921 : speckleLifetimes( mnMax + m, mnMax + n ) = tau;
1922 : speckleLifetimes( mnMax - m, mnMax - n ) = tau;
1923 : }
1924 :
1925 : if( doLP )
1926 : {
1927 : if( gopt_lp > 0 )
1928 : {
1929 : for( size_t i = 0; i < tfreq.size(); ++i )
1930 : {
1931 : ETFxn[i] = go_lp.clETF( i, gopt_lp );
1932 : NTFxn[i] = go_lp.clNTF( i, gopt_lp );
1933 : }
1934 : }
1935 : else
1936 : {
1937 : for( size_t i = 0; i < tfreq.size(); ++i )
1938 : {
1939 : ETFxn[i] = 1;
1940 : NTFxn[i] = 0;
1941 : }
1942 : }
1943 :
1944 : if( writeXfer )
1945 : {
1946 : std::string tfOutFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] );
1947 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
1948 :
1949 : std::string etfOutFile = std::format( "{}/etf_{}_{}.binv", tfOutFile, m, n );
1950 : // tfOutFile + "etf_" + ioutils::convert ToString( m ) + '_' +
1951 : // ioutils::convert ToString( n ) + ".binv";
1952 : ioutils::writeBinVector( etfOutFile, ETFxn );
1953 :
1954 : std::string ntfOutFile = std::format( "{}/ntf_{}_{}.binv", tfOutFile, m, n );
1955 : // tfOutFile + "ntf_" + ioutils::convert ToString( m ) + '_' +
1956 : // ioutils::convert ToString( n ) + ".binv";
1957 : ioutils::writeBinVector( ntfOutFile, NTFxn );
1958 :
1959 : if( i == 0 ) // Write freq on the first one
1960 : {
1961 : std::string fOutFile = tfOutFile + "freq.binv";
1962 : ioutils::writeBinVector( fOutFile, tfreq );
1963 : }
1964 : }
1965 :
1966 : if( lifetimeTrials > 0 )
1967 : {
1968 : speckleAmpPSD( spfreq, sppsd, tfreq, tPSDp, ETFxn, tPSDn, NTFxn, lifetimeTrials );
1969 : realT spvar = sigproc::psdVar( spfreq, sppsd );
1970 :
1971 : realT splifeT = 100.0;
1972 : realT error;
1973 :
1974 : realT tau = pvm( error, spfreq, sppsd, splifeT ) * ( splifeT ) / spvar;
1975 :
1976 : speckleLifetimes_lp( mnMax + m, mnMax + n ) = tau;
1977 : speckleLifetimes_lp( mnMax - m, mnMax - n ) = tau;
1978 : }
1979 : } // if(doLP)
1980 : } // if( (lifetimeTrials > 0 || writeXfer) && ( uncontrolledLifetimes || inside ))
1981 : //**>
1982 :
1983 : // Calculate the controlled PSDs and output
1984 : if( writePSDs )
1985 : {
1986 : std::string psdOutFile =
1987 : std::format( "{}/outputPSDs_{}_si/psd_{}_{}.binv", dir, mags[s], m, n );
1988 :
1989 : // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/";
1990 : // psdOutFile +=
1991 : // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
1992 :
1993 : std::vector<realT> psdOut( tPSDp.size() + tPSDpHF.size() );
1994 :
1995 : // Calculate the output PSD if gains are applied
1996 : if( gopt > 0 )
1997 : {
1998 : realT ETF, NTF;
1999 :
2000 : for( size_t i = 0; i < tfreq.size(); ++i )
2001 : {
2002 : go_si.clTF2( ETF, NTF, i, gopt );
2003 : psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
2004 : }
2005 :
2006 : for( size_t i = 0; i < tPSDpHF.size(); ++i )
2007 : {
2008 : psdOut[tfreq.size() + i] = tPSDpHF[i];
2009 : }
2010 : }
2011 : else // otherwise just copy
2012 : {
2013 : for( size_t i = 0; i < tfreq.size(); ++i )
2014 : {
2015 : psdOut[i] = tPSDp[i];
2016 : }
2017 :
2018 : for( size_t i = 0; i < tPSDpHF.size(); ++i )
2019 : {
2020 : psdOut[tfreq.size() + i] = tPSDpHF[i];
2021 : }
2022 : }
2023 :
2024 : ioutils::writeBinVector( psdOutFile, psdOut );
2025 :
2026 : if( i == 0 ) // Write freq on the first one
2027 : {
2028 : psdOutFile = std::format( "{}/outputPSDs_{}_si/freq.binv", dir, mags[s] );
2029 : // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_si/freq.binv";
2030 : // ioutils::writeBinVector( psdOutFile, tfreqHF );
2031 : }
2032 :
2033 : if( doLP )
2034 : {
2035 : std::string psdOutFile =
2036 : std::format( "{}/outputPSDs_{}_lp/psd_{}_{}.binv", dir, mags[s], m, n );
2037 :
2038 : // psdOutFile = dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) + "_lp/";
2039 : // psdOutFile +=
2040 : // "psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) +
2041 : // ".binv";
2042 :
2043 : // Calculate the output PSD if gains are applied
2044 : if( gopt_lp > 0 )
2045 : {
2046 : realT ETF, NTF;
2047 :
2048 : for( size_t i = 0; i < tfreq.size(); ++i )
2049 : {
2050 : go_lp.clTF2( ETF, NTF, i, gopt_lp );
2051 : psdOut[i] = tPSDp[i] * ETF + tPSDn[i] * NTF;
2052 : }
2053 : for( size_t i = 0; i < tPSDpHF.size(); ++i )
2054 : {
2055 : psdOut[tfreq.size() + i] = tPSDpHF[i];
2056 : }
2057 : }
2058 : else // otherwise just copy
2059 : {
2060 : for( size_t i = 0; i < tfreq.size(); ++i )
2061 : {
2062 : psdOut[i] = tPSDp[i];
2063 : }
2064 :
2065 : for( size_t i = 0; i < tPSDpHF.size(); ++i )
2066 : {
2067 : psdOut[tfreq.size() + i] = tPSDpHF[i];
2068 : }
2069 : }
2070 :
2071 : ioutils::writeBinVector( psdOutFile, psdOut );
2072 :
2073 : if( i == 0 )
2074 : {
2075 : psdOutFile = std::format( "{}/outputPSDs_{}_lp/freq.binv", dir, mags[s] );
2076 : // psdOutFile =
2077 : // dir + "/" + "outputPSDs_" + ioutils::convert ToString( mags[s] ) +
2078 : // "_lp/freq.binv";
2079 : // ioutils::writeBinVector( psdOutFile, tfreq );
2080 : }
2081 : }
2082 : }
2083 : }
2084 : watcher.incrementAndOutputStatus();
2085 :
2086 : } // omp for i..nModes
2087 : } // omp Parallel
2088 :
2089 0 : if( analysisStatus.load() != static_cast<int>( error_t::noerror ) )
2090 : {
2091 0 : return analysisStatus.load();
2092 : }
2093 :
2094 0 : Eigen::Array<realT, -1, -1> cim;
2095 :
2096 0 : fits::fitsFile<realT> ff;
2097 0 : std::string fn = std::format( "{}/gainmap_{}_si.fits", dir, mags[s] );
2098 : // dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2099 0 : ff.write( fn, gains );
2100 :
2101 0 : fn = std::format( "{}/varmap_{}_si.fits", dir, mags[s] );
2102 : // dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2103 0 : ff.write( fn, vars );
2104 :
2105 0 : cim = vars;
2106 :
2107 0 : realT Ssi = exp( -1 * cim.sum() );
2108 0 : S_si.push_back( strehl );
2109 0 : cim /= strehl;
2110 :
2111 0 : fn = std::format( "{}/contrast_{}_si.fits", dir, mags[s] );
2112 : // dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2113 0 : ff.write( fn, cim );
2114 :
2115 0 : if( lifetimeTrials > 0 )
2116 : {
2117 0 : fn = std::format( "{}/speckleLifetimes_{}_si.fits", dir, mags[s] );
2118 : // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2119 0 : ff.write( fn, speckleLifetimes );
2120 : }
2121 :
2122 0 : if( doLP )
2123 : {
2124 0 : fn = std::format( "{}/gainmap_{}_lp.fits", dir, mags[s] );
2125 : // dir + "/gainmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2126 0 : ff.write( fn, gains_lp );
2127 :
2128 0 : fn = std::format( "{}/lpcmap_{}_lp.fits", dir, mags[s] );
2129 : // dir + "/lpcmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2130 0 : ff.write( fn, lpC );
2131 :
2132 0 : fn = std::format( "{}/varmap_{}_lp.fits", dir, mags[s] );
2133 : // dir + "/varmap_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2134 0 : ff.write( fn, vars_lp );
2135 :
2136 0 : cim = vars_lp;
2137 :
2138 : // Scale Strehl by the ratio of total variance
2139 0 : realT Slp = strehl * exp( -1 * cim.sum() ) /
2140 : Ssi; // This is a hack until we do a real fitting error calculation or something
2141 0 : S_lp.push_back( Slp );
2142 0 : cim /= Slp;
2143 :
2144 0 : fn = std::format( "{}/contrast_{}_lp.fits", dir, mags[s] );
2145 : // dir + "/contrast_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2146 0 : ff.write( fn, cim );
2147 :
2148 0 : if( lifetimeTrials > 0 )
2149 : {
2150 0 : fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] );
2151 : // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2152 0 : ff.write( fn, speckleLifetimes_lp );
2153 : }
2154 : }
2155 :
2156 : } // s (mag)
2157 :
2158 0 : fn = dir + "/strehl_si.txt";
2159 0 : fout.open( fn );
2160 0 : for( size_t i = 0; i < mags.size(); ++i )
2161 : {
2162 0 : fout << mags[i] << " " << S_si[i] << "\n";
2163 : }
2164 :
2165 0 : fout.close();
2166 :
2167 0 : if( doLP )
2168 : {
2169 0 : fn = dir + "/strehl_lp.txt";
2170 0 : fout.open( fn );
2171 0 : for( size_t i = 0; i < mags.size(); ++i )
2172 : {
2173 0 : fout << mags[i] << " " << S_lp[i] << "\n";
2174 : }
2175 :
2176 0 : fout.close();
2177 : }
2178 :
2179 0 : return 0;
2180 0 : }
2181 :
2182 : template <typename realT, typename aosysT>
2183 0 : int fourierTemporalPSD<realT, aosysT>::intensityPSD(
2184 : const std::string &subDir, // sub-directory of psdDir which contains the controlled system results,
2185 : // and where the lifetimes will be written.
2186 : const std::string &psdDir, // directory containing the PSDS
2187 : const std::string &CvdPath, // path to the covariance decomposition
2188 : int mnMax,
2189 : int mnCon,
2190 : std::vector<realT> &mags,
2191 : int lifetimeTrials,
2192 : bool writePSDs )
2193 : {
2194 :
2195 0 : std::string dir = psdDir + "/" + subDir;
2196 :
2197 : /*** Dump Params to file ***/
2198 0 : mkdir( dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
2199 :
2200 0 : std::ofstream fout;
2201 0 : std::string fn = dir + '/' + "splife_params.txt";
2202 0 : fout.open( fn );
2203 :
2204 0 : fout << "#---------------------------\n";
2205 0 : m_aosys->dumpAOSystem( fout );
2206 0 : fout << "#---------------------------\n";
2207 0 : fout << "# Analysis Parameters\n";
2208 0 : fout << "# mnMax = " << mnMax << '\n';
2209 0 : fout << "# mnCon = " << mnCon << '\n';
2210 0 : fout << "# mags = ";
2211 0 : for( size_t i = 0; i < mags.size() - 1; ++i )
2212 0 : fout << mags[i] << ",";
2213 0 : fout << mags[mags.size() - 1] << '\n';
2214 0 : fout << "# lifetimeTrials = " << lifetimeTrials << '\n';
2215 : // fout << "# uncontrolledLifetimes = " << uncontrolledLifetimes << '\n';
2216 0 : fout << "# writePSDs = " << std::boolalpha << writePSDs << '\n';
2217 :
2218 0 : fout.close();
2219 :
2220 0 : realT fs = 1.0 / m_aosys->tauWFS();
2221 0 : realT tauWFS = m_aosys->tauWFS();
2222 0 : realT deltaTau = m_aosys->deltaTau();
2223 :
2224 : //** Get the Fourier Grid
2225 0 : std::vector<sigproc::fourierModeDef> fms;
2226 :
2227 0 : sigproc::makeFourierModeFreqs_Rect( fms, 2 * mnMax );
2228 0 : size_t nModes = 0.5 * fms.size();
2229 0 : std::cerr << "nModes: " << nModes << " (" << fms.size() << ")\n";
2230 :
2231 0 : Eigen::Array<realT, -1, -1> speckleLifetimes;
2232 0 : Eigen::Array<realT, -1, -1> speckleLifetimes_lp;
2233 :
2234 0 : speckleLifetimes.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2235 0 : speckleLifetimes( mnMax, mnMax ) = 0;
2236 :
2237 0 : speckleLifetimes_lp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2238 0 : speckleLifetimes_lp( mnMax, mnMax ) = 0;
2239 :
2240 : /*********************************************************************/
2241 : // 0) Get the frequency grid, and nyquist limit it to f_s/2
2242 : /*********************************************************************/
2243 :
2244 0 : std::vector<realT> tfreq;
2245 0 : std::vector<realT> tPSDp; // The open-loop OPD PSD
2246 0 : std::vector<realT> tPSDn; // The open-loop WFS noise PSD
2247 0 : std::vector<complexT> tETF;
2248 0 : std::vector<complexT> tNTF;
2249 :
2250 0 : if( getGridFreq( tfreq, psdDir ) < 0 )
2251 0 : return -1;
2252 :
2253 0 : size_t imax = 0;
2254 0 : while( tfreq[imax] <= 0.5 * fs )
2255 : {
2256 0 : ++imax;
2257 0 : if( imax > tfreq.size() - 1 )
2258 0 : break;
2259 : }
2260 :
2261 0 : if( imax < tfreq.size() - 1 && tfreq[imax] <= 0.5 * fs * ( 1.0 + 1e-7 ) )
2262 0 : ++imax;
2263 :
2264 0 : tfreq.erase( tfreq.begin() + imax, tfreq.end() );
2265 :
2266 : // Now allocate memory
2267 0 : tPSDn.resize( tfreq.size() );
2268 0 : std::vector<std::vector<realT>> sqrtOPDPSD;
2269 0 : sqrtOPDPSD.resize( nModes );
2270 :
2271 0 : std::vector<std::vector<realT>> opdPSD;
2272 0 : opdPSD.resize( nModes );
2273 :
2274 0 : std::vector<realT> psd2sided;
2275 :
2276 : // Store the mode variance for later normalization
2277 0 : std::vector<realT> modeVar;
2278 0 : modeVar.resize( nModes );
2279 :
2280 : /*********************************************************************/
2281 : // 1) Read in each PSD, and load it into the array in FFT order
2282 : /*********************************************************************/
2283 :
2284 0 : for( size_t i = 0; i < nModes; ++i )
2285 : {
2286 : // Determine the spatial frequency at this step
2287 0 : int m = fms[2 * i].m;
2288 0 : int n = fms[2 * i].n;
2289 :
2290 : //**< Get the open-loop turb. PSD
2291 0 : if( getGridPSD( tPSDp, psdDir, m, n ) < 0 )
2292 0 : return -1;
2293 0 : tPSDp.erase( tPSDp.begin() + imax, tPSDp.end() ); // Nyquist limit
2294 0 : modeVar[i] = sigproc::psdVar( tfreq, tPSDp );
2295 :
2296 : // And now normalize
2297 0 : sigproc::normPSD( tPSDp, tfreq, 1.0 ); // Normalize
2298 0 : sigproc::augment1SidedPSD( psd2sided, tPSDp, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
2299 :
2300 0 : opdPSD[i].resize( psd2sided.size() );
2301 0 : sqrtOPDPSD[i].resize( psd2sided.size() );
2302 :
2303 0 : for( size_t j = 0; j < psd2sided.size(); ++j )
2304 : {
2305 0 : opdPSD[i][j] = psd2sided[j] * modeVar[i];
2306 0 : sqrtOPDPSD[i][j] = sqrt( psd2sided[j] ); // Store the square-root for later
2307 : }
2308 : }
2309 :
2310 0 : size_t sz2Sided = psd2sided.size();
2311 :
2312 0 : std::vector<realT> freq2sided;
2313 0 : freq2sided.resize( sz2Sided );
2314 0 : sigproc::augment1SidedPSDFreq( freq2sided, tfreq );
2315 :
2316 0 : tPSDp.resize( tfreq.size() );
2317 0 : tETF.resize( tfreq.size() );
2318 0 : tNTF.resize( tfreq.size() );
2319 :
2320 0 : std::vector<std::vector<realT>> sqrtNPSD;
2321 0 : sqrtNPSD.resize( nModes );
2322 :
2323 0 : std::vector<realT> noiseVar;
2324 0 : noiseVar.resize( nModes );
2325 :
2326 0 : std::vector<std::vector<complexT>> ETFsi;
2327 0 : std::vector<std::vector<complexT>> ETFlp;
2328 0 : ETFsi.resize( nModes );
2329 0 : ETFlp.resize( nModes );
2330 :
2331 0 : std::vector<std::vector<complexT>> NTFsi;
2332 0 : std::vector<std::vector<complexT>> NTFlp;
2333 0 : NTFsi.resize( nModes );
2334 0 : NTFlp.resize( nModes );
2335 :
2336 0 : std::string tfInFile;
2337 0 : std::string etfInFile;
2338 0 : std::string ntfInFile;
2339 :
2340 0 : improc::eigenImage<realT> Cvd;
2341 0 : fits::fitsFile<realT> ff;
2342 0 : ff.read( Cvd, CvdPath );
2343 :
2344 0 : std::vector<std::complex<realT>> tPSDc, psd2sidedc;
2345 :
2346 : /*********************************************************************/
2347 : // 2) Analyze each star magnitude
2348 : /*********************************************************************/
2349 0 : ipc::ompLoopWatcher<> watcher( lifetimeTrials * mags.size(), std::cout );
2350 0 : for( size_t s = 0; s < mags.size(); ++s )
2351 : {
2352 : /*********************************************************************/
2353 : // 2.0) Read in the transfer functions for each mode
2354 : /*********************************************************************/
2355 0 : for( size_t i = 0; i < nModes; ++i )
2356 : {
2357 : // Determine the spatial frequency at this step
2358 0 : int m = fms[2 * i].m;
2359 0 : int n = fms[2 * i].n;
2360 :
2361 : //**< Determine if we're inside the hardwarecontrol limit
2362 0 : bool inside = false;
2363 :
2364 0 : if( m_aosys->circularLimit() )
2365 : {
2366 0 : if( m * m + n * n <= mnCon * mnCon )
2367 0 : inside = true;
2368 : }
2369 : else
2370 : {
2371 0 : if( fabs( m ) <= mnCon && fabs( n ) <= mnCon )
2372 0 : inside = true;
2373 : }
2374 :
2375 : // Get the WFS noise PSD (which is already resized to match tfreq)
2376 0 : wfsNoisePSD<realT>( tPSDn,
2377 0 : m_aosys->beta_p( m, n ),
2378 0 : m_aosys->Fg( mags[s] ),
2379 : tauWFS,
2380 0 : m_aosys->npix_wfs( (size_t)0 ),
2381 0 : m_aosys->Fbg( (size_t)0 ),
2382 0 : m_aosys->ron_wfs( (size_t)0 ) );
2383 0 : sigproc::augment1SidedPSD( psd2sided, tPSDn, !( tfreq[0] == 0 ) ); // Convert to FFT storage order
2384 :
2385 : // Pre-calculate the variance of the noise for later use
2386 0 : noiseVar[i] = sigproc::psdVar( tfreq, tPSDn );
2387 :
2388 0 : sqrtNPSD[i].resize( psd2sided.size() );
2389 0 : for( size_t j = 0; j < psd2sided.size(); ++j )
2390 0 : sqrtNPSD[i][j] = sqrt( psd2sided[j] );
2391 :
2392 0 : ETFsi[i].resize( sz2Sided );
2393 0 : ETFlp[i].resize( sz2Sided );
2394 0 : NTFsi[i].resize( sz2Sided );
2395 0 : NTFlp[i].resize( sz2Sided );
2396 :
2397 0 : if( inside )
2398 : {
2399 0 : tfInFile = std::format( "{}/outputTF_{}_si", dir, mags[s] );
2400 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_si/";
2401 :
2402 0 : etfInFile = std::format( "{}/etf_{}_{}.binv", tfInFile, m, n );
2403 : // tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2404 0 : ioutils::readBinVector( tPSDc, etfInFile );
2405 0 : sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
2406 0 : for( size_t j = 0; j < psd2sidedc.size(); ++j )
2407 0 : ETFsi[i][j] = psd2sidedc[j];
2408 :
2409 0 : ntfInFile = std::format( "{}/ntf_{}_{}.binv", tfInFile, m, n );
2410 : // tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2411 0 : ioutils::readBinVector( tPSDc, ntfInFile );
2412 0 : sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
2413 0 : for( size_t j = 0; j < psd2sidedc.size(); ++j )
2414 0 : NTFsi[i][j] = psd2sidedc[j];
2415 :
2416 0 : tfInFile = std::format( "{}/outputTF_{}_lp", dir, mags[s] );
2417 : // dir + "/" + "outputTF_" + ioutils::convert ToString( mags[s] ) + "_lp/";
2418 :
2419 0 : etfInFile = std::format( "{}/etf_{}_{}.binv", tfInFile, m, n );
2420 : // tfInFile + "etf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2421 0 : ioutils::readBinVector( tPSDc, etfInFile );
2422 0 : sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
2423 0 : for( size_t j = 0; j < psd2sidedc.size(); ++j )
2424 0 : ETFlp[i][j] = psd2sidedc[j];
2425 :
2426 0 : ntfInFile = std::format( "{}/ntf_{}_{}.binv", tfInFile, m, n );
2427 : // tfInFile + "ntf_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2428 0 : ioutils::readBinVector( tPSDc, ntfInFile );
2429 0 : sigproc::augment1SidedPSD( psd2sidedc, tPSDc, !( tfreq[0] == 0 ), 1 ); // Convert to FFT storage order
2430 0 : for( size_t j = 0; j < psd2sidedc.size(); ++j )
2431 0 : NTFlp[i][j] = psd2sidedc[j];
2432 : }
2433 : else
2434 : {
2435 0 : for( int q = 0; q < ETFsi.size(); ++q )
2436 : {
2437 0 : ETFsi[i][q] = 1;
2438 0 : NTFsi[i][q] = 0;
2439 0 : ETFlp[i][q] = 1;
2440 0 : NTFlp[i][q] = 0;
2441 : }
2442 : }
2443 : }
2444 :
2445 0 : sigproc::averagePeriodogram<realT> tavgPgram(
2446 : sz2Sided / 1.,
2447 0 : 1 / fs ); // this is just to get the size right, per-thread instances below
2448 0 : std::vector<std::vector<realT>> spPSDs;
2449 0 : spPSDs.resize( nModes );
2450 0 : for( size_t pp = 0; pp < spPSDs.size(); ++pp )
2451 : {
2452 0 : spPSDs[pp].resize( tavgPgram.size() );
2453 0 : for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2454 0 : spPSDs[pp][nn] = 0;
2455 : }
2456 :
2457 0 : std::vector<std::vector<realT>> spPSDslp;
2458 0 : spPSDslp.resize( nModes );
2459 0 : for( size_t pp = 0; pp < spPSDslp.size(); ++pp )
2460 : {
2461 0 : spPSDslp[pp].resize( tavgPgram.size() );
2462 0 : for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2463 0 : spPSDslp[pp][nn] = 0;
2464 : }
2465 :
2466 0 : #pragma omp parallel
2467 : {
2468 : // Normally distributed random numbers
2469 : math::normDistT<realT> normVar;
2470 : normVar.seed();
2471 :
2472 : // FFTs for going to Fourier domain and back to time domain.
2473 : math::ft::fftT<realT, std::complex<realT>, 1, 0> fftF( sqrtOPDPSD[0].size() );
2474 : math::ft::fftT<std::complex<realT>, realT, 1, 0> fftB( sqrtOPDPSD[0].size(), math::ft::dir::backward );
2475 :
2476 : // Working memory
2477 : std::vector<std::complex<realT>> tform1( sqrtOPDPSD[0].size() );
2478 : std::vector<std::complex<realT>> tform2( sqrtOPDPSD[0].size() );
2479 : std::vector<std::complex<realT>> Ntform1( sqrtOPDPSD[0].size() );
2480 : std::vector<std::complex<realT>> Ntform2( sqrtOPDPSD[0].size() );
2481 :
2482 : std::vector<std::complex<realT>> tform1lp( sqrtOPDPSD[0].size() );
2483 : std::vector<std::complex<realT>> tform2lp( sqrtOPDPSD[0].size() );
2484 : std::vector<std::complex<realT>> Ntform1lp( sqrtOPDPSD[0].size() );
2485 : std::vector<std::complex<realT>> Ntform2lp( sqrtOPDPSD[0].size() );
2486 :
2487 : // OPD-PSD filter
2488 : sigproc::psdFilter<realT, 1> pfilt;
2489 : pfilt.psdSqrt( &sqrtOPDPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
2490 :
2491 : // Noise-PSD filter
2492 : sigproc::psdFilter<realT, 1> nfilt;
2493 : nfilt.psdSqrt( &sqrtNPSD[0], tfreq[1] - tfreq[0] ); // Pre-configure
2494 :
2495 : // The h time-series
2496 : std::vector<std::vector<realT>> hts;
2497 : hts.resize( 2 * nModes );
2498 :
2499 : // The correlated h time-series
2500 : std::vector<std::vector<realT>> htsCorr;
2501 : htsCorr.resize( 2 * nModes );
2502 :
2503 : for( size_t pp = 0; pp < hts.size(); ++pp )
2504 : {
2505 : hts[pp].resize( sqrtOPDPSD[0].size() );
2506 : htsCorr[pp].resize( sqrtOPDPSD[0].size() );
2507 : }
2508 :
2509 : // The noise time-serieses
2510 : std::vector<realT> N_n;
2511 : N_n.resize( sz2Sided );
2512 :
2513 : std::vector<realT> N_nm;
2514 : N_nm.resize( sz2Sided );
2515 :
2516 : // Periodogram averager
2517 : sigproc::averagePeriodogram<realT> avgPgram( sz2Sided / 1., 1 / fs ); //, 0, 1);
2518 : avgPgram.win( sigproc::window::hann );
2519 :
2520 : // The periodogram output
2521 : std::vector<realT> tpgram( avgPgram.size() );
2522 :
2523 : // Holds the speckle time-series
2524 : improc::eigenCube<realT> spTS;
2525 : spTS.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
2526 :
2527 : improc::eigenCube<realT> spTSlp;
2528 : spTSlp.resize( 2 * mnMax + 1, 2 * mnMax + 1, tform1.size() );
2529 :
2530 : // Here's where the big loop of n-trials should start
2531 : #pragma omp for
2532 : for( int zz = 0; zz < lifetimeTrials; ++zz )
2533 : {
2534 : std::complex<realT> scale = exp( std::complex<realT>( 0, math::half_pi<realT>() ) ) /
2535 : std::complex<realT>( ( tform1.size() ), 0 );
2536 :
2537 : /*********************************************************************/
2538 : // 2.1) Generate filtered noise for each mode, with temporal phase shifts at each spatial frequency
2539 : /*********************************************************************/
2540 : for( size_t pp = 0; pp < nModes; ++pp )
2541 : {
2542 : // Fill in standard normal noise
2543 : for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
2544 : {
2545 : hts[2 * pp][nn] = normVar;
2546 : }
2547 :
2548 : // Set sqrt(PSD), just a pointer switch
2549 : pfilt.psdSqrt( &sqrtOPDPSD[pp], tfreq[1] - tfreq[0] );
2550 :
2551 : // And now filter the noise to a time-series of h
2552 : pfilt( hts[2 * pp] );
2553 :
2554 : /**/
2555 : // Then construct 2nd mode with temporal shift
2556 : fftF( tform1.data(), hts[2 * pp].data() );
2557 :
2558 : // Apply the phase shift to form the 2nd time series
2559 : for( size_t nn = 0; nn < hts[2 * pp].size(); ++nn )
2560 : tform1[nn] = tform1[nn] * scale;
2561 :
2562 : fftB( hts[2 * pp + 1].data(), tform1.data() );
2563 : /**/
2564 : }
2565 : //** At this point we have correlated time-series, with the correct temporal PSD, but not yet spatially
2566 : // correlated
2567 :
2568 : /*********************************************************************/
2569 : // 2.2) Correlate the time-series for each mode
2570 : /*********************************************************************/
2571 : // #pragma omp parallel for
2572 : for( size_t pp = 0; pp < hts.size(); ++pp )
2573 : {
2574 : for( size_t nn = 0; nn < hts[0].size(); ++nn )
2575 : {
2576 : htsCorr[pp][nn] = 0;
2577 : }
2578 :
2579 : for( size_t qq = 0; qq <= pp; ++qq )
2580 : {
2581 : realT cvd = Cvd( qq, pp );
2582 : realT *d1 = htsCorr[pp].data();
2583 : realT *d2 = hts[qq].data();
2584 : for( size_t nn = 0; nn < hts[0].size(); ++nn )
2585 : {
2586 : d1[nn] += d2[nn] * cvd;
2587 : }
2588 : }
2589 : }
2590 :
2591 : /*
2592 : for(size_t pp=0; pp < hts.size(); ++pp)
2593 : {
2594 : for(size_t nn=0; nn< hts[0].size(); ++nn)
2595 : {
2596 : htsCorr[pp][nn] = hts[pp][nn];
2597 : }
2598 : }*/
2599 :
2600 : /*********************************************************************/
2601 : // 2.2.a) Re-normalize b/c the correlation step above does not result in correct variances
2602 : ///\todo should be able to scale the covar by r0, and possibly D
2603 : /*********************************************************************/
2604 : for( size_t pp = 0; pp < nModes; ++pp )
2605 : {
2606 : math::vectorMeanSub( htsCorr[2 * pp] );
2607 : math::vectorMeanSub( htsCorr[2 * pp + 1] );
2608 :
2609 : realT var = math::vectorVariance( htsCorr[2 * pp] );
2610 : realT norm = sqrt( modeVar[pp] / var );
2611 : for( size_t nn = 0; nn < htsCorr[2 * pp].size(); ++nn )
2612 : htsCorr[2 * pp][nn] *= norm;
2613 :
2614 : var = math::vectorVariance( htsCorr[2 * pp + 1] );
2615 : norm = sqrt( modeVar[pp] / var );
2616 : for( size_t nn = 0; nn < htsCorr[2 * pp + 1].size(); ++nn )
2617 : htsCorr[2 * pp + 1][nn] *= norm;
2618 : }
2619 :
2620 : scale = std::complex<realT>( tform1.size(), 0 );
2621 :
2622 : /*********************************************************************/
2623 : // 2.3) Generate speckle intensity time-series
2624 : /*********************************************************************/
2625 : for( size_t pp = 0; pp < nModes; ++pp )
2626 : {
2627 : // Now we take them back to the FD and apply the xfer
2628 : // and add the noise
2629 :
2630 : fftF( tform1.data(), htsCorr[2 * pp].data() );
2631 : fftF( tform2.data(), htsCorr[2 * pp + 1].data() );
2632 :
2633 : // Make some noise
2634 : for( int nn = 0; nn < sz2Sided; ++nn )
2635 : {
2636 : N_n[nn] = normVar;
2637 : N_nm[nn] = normVar;
2638 : }
2639 :
2640 : // Filter it
2641 : // Set sqrt(PSD), just a pointer switch
2642 : pfilt.psdSqrt( &sqrtNPSD[pp], tfreq[1] - tfreq[0] );
2643 : nfilt.filter( N_n );
2644 : nfilt.filter( N_nm );
2645 :
2646 : // Normalize it
2647 : realT Nactvar = 0.5 * ( math::vectorVariance( N_n ) + math::vectorVariance( N_nm ) );
2648 : realT norm = sqrt( noiseVar[pp] / Nactvar );
2649 : for( size_t q = 0; q < N_n.size(); ++q )
2650 : N_n[q] *= norm;
2651 : for( size_t q = 0; q < N_nm.size(); ++q )
2652 : N_nm[q] *= norm;
2653 :
2654 : // And move them to the Fourier domain
2655 : fftF( Ntform1.data(), N_n.data() );
2656 : fftF( Ntform2.data(), N_nm.data() );
2657 :
2658 : // Apply the closed loop transfers
2659 : for( size_t mm = 0; mm < tform1.size(); ++mm )
2660 : {
2661 : // Apply the augmented ETF to two time-series
2662 : tform1lp[mm] = tform1[mm] * ETFlp[pp][mm] / scale;
2663 : tform2lp[mm] = tform2[mm] * ETFlp[pp][mm] / scale;
2664 :
2665 : Ntform1lp[mm] = Ntform1[mm] * NTFlp[pp][mm] / scale;
2666 : Ntform2lp[mm] = Ntform2[mm] * NTFlp[pp][mm] / scale;
2667 :
2668 : tform1[mm] *= ETFsi[pp][mm] / scale;
2669 : tform2[mm] *= ETFsi[pp][mm] / scale;
2670 :
2671 : Ntform1[mm] *= NTFsi[pp][mm] / scale;
2672 : Ntform2[mm] *= NTFsi[pp][mm] / scale;
2673 : }
2674 :
2675 : // And make the speckle TS
2676 : int m = fms[2 * pp].m;
2677 : int n = fms[2 * pp].n;
2678 :
2679 : //<<<<<<<<****** Transform back to the time domain.
2680 : fftB( htsCorr[2 * pp].data(), tform1.data() );
2681 : fftB( htsCorr[2 * pp + 1].data(), tform2.data() );
2682 : fftB( N_n.data(), Ntform1.data() );
2683 : fftB( N_nm.data(), Ntform2.data() );
2684 :
2685 : for( int i = 0; i < hts[2 * pp].size(); ++i )
2686 : {
2687 : realT h1 = htsCorr[2 * pp][i] + N_n[i];
2688 : realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2689 :
2690 : spTS.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2691 : spTS.image( i )( mnMax - m, mnMax - n ) = spTS.image( i )( mnMax + m, mnMax + n );
2692 : }
2693 :
2694 : fftB( htsCorr[2 * pp].data(), tform1lp.data() );
2695 : fftB( htsCorr[2 * pp + 1].data(), tform2lp.data() );
2696 : fftB( N_n.data(), Ntform1lp.data() );
2697 : fftB( N_nm.data(), Ntform2lp.data() );
2698 :
2699 : for( int i = 0; i < hts[2 * pp].size(); ++i )
2700 : {
2701 : realT h1 = htsCorr[2 * pp][i] + N_n[i];
2702 : realT h2 = htsCorr[2 * pp + 1][i] + N_nm[i];
2703 :
2704 : spTSlp.image( i )( mnMax + m, mnMax + n ) = ( pow( h1, 2 ) + pow( h2, 2 ) );
2705 : spTSlp.image( i )( mnMax - m, mnMax - n ) = spTSlp.image( i )( mnMax + m, mnMax + n );
2706 : }
2707 : }
2708 :
2709 : /*********************************************************************/
2710 : // 2.5) Calculate speckle PSD for each mode
2711 : /*********************************************************************/
2712 : std::vector<realT> speckAmp( spTS.planes() );
2713 : std::vector<realT> speckAmplp( spTS.planes() );
2714 :
2715 : for( size_t pp = 0; pp < nModes; ++pp )
2716 : {
2717 : int m = fms[2 * pp].m;
2718 : int n = fms[2 * pp].n;
2719 :
2720 : realT mn = 0;
2721 : realT mnlp = 0;
2722 : for( int i = 0; i < spTS.planes(); ++i )
2723 : {
2724 : speckAmp[i] = spTS.image( i )( mnMax + m, mnMax + n );
2725 : speckAmplp[i] = spTSlp.image( i )( mnMax + m, mnMax + n );
2726 :
2727 : mn += speckAmp[i];
2728 : mnlp += speckAmplp[i];
2729 : }
2730 : mn /= speckAmp.size();
2731 : mnlp /= speckAmplp.size();
2732 :
2733 : // mean subtract
2734 : for( int i = 0; i < speckAmp.size(); ++i )
2735 : speckAmp[i] -= mn;
2736 : for( int i = 0; i < speckAmplp.size(); ++i )
2737 : speckAmplp[i] -= mnlp;
2738 :
2739 : // Calculate PSD of the speckle amplitude
2740 : avgPgram( tpgram, speckAmp );
2741 : for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2742 : spPSDs[pp][nn] += tpgram[nn];
2743 :
2744 : avgPgram( tpgram, speckAmplp );
2745 : for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2746 : spPSDslp[pp][nn] += tpgram[nn];
2747 : }
2748 :
2749 : watcher.incrementAndOutputStatus();
2750 : } // for(int zz=0; zz<lifetimeTrials; ++zz)
2751 : } // #pragma omp parallel
2752 :
2753 0 : std::vector<realT> spFreq( spPSDs[0].size() );
2754 0 : for( size_t nn = 0; nn < spFreq.size(); ++nn )
2755 0 : spFreq[nn] = tavgPgram[nn];
2756 :
2757 0 : improc::eigenImage<realT> taus, tauslp;
2758 0 : taus.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2759 0 : tauslp.resize( 2 * mnMax + 1, 2 * mnMax + 1 );
2760 :
2761 0 : improc::eigenCube<realT> imc, imclp;
2762 0 : std::vector<realT> clPSD;
2763 :
2764 0 : if( writePSDs )
2765 : {
2766 0 : imc.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2767 0 : imclp.resize( 2 * mnMax + 1, 2 * mnMax + 1, spPSDs[0].size() );
2768 0 : clPSD.resize( sz2Sided );
2769 : }
2770 :
2771 0 : sigproc::psdVarMean<sigproc::psdVarMeanParams<realT>> pvm;
2772 : /*********************************************************************/
2773 : // 3.0) Calculate lifetimes from the PSDs
2774 : /*********************************************************************/
2775 0 : for( size_t pp = 0; pp < nModes; ++pp )
2776 : {
2777 0 : spPSDs[pp][0] = spPSDs[pp][1]; // deal with under-estimated mean.
2778 0 : spPSDslp[pp][0] = spPSDslp[pp][1]; // deal with under-estimated mean.
2779 :
2780 0 : int m = fms[2 * pp].m;
2781 0 : int n = fms[2 * pp].n;
2782 :
2783 : realT var;
2784 0 : if( writePSDs ) // Have to normalize the intensity for some reason if we want to use the PSDs
2785 : {
2786 0 : for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2787 : {
2788 0 : spPSDs[pp][nn] /= lifetimeTrials;
2789 : }
2790 :
2791 0 : for( size_t nn = 0; nn < sz2Sided; ++nn )
2792 : {
2793 0 : clPSD[nn] =
2794 0 : opdPSD[pp][nn] * norm( ETFsi[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFsi[pp][nn] );
2795 : }
2796 :
2797 0 : var = sigproc::psdVar( freq2sided, clPSD );
2798 :
2799 0 : realT pvar = sigproc::psdVar( spFreq, spPSDs[pp] );
2800 :
2801 0 : for( size_t nn = 0; nn < spPSDs[pp].size(); ++nn )
2802 : {
2803 0 : spPSDs[pp][nn] *= var / pvar;
2804 0 : imc.image( nn )( mnMax + m, mnMax + n ) = spPSDs[pp][nn];
2805 0 : imc.image( nn )( mnMax - m, mnMax - n ) = spPSDs[pp][nn];
2806 : }
2807 :
2808 : // lp
2809 0 : for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2810 : {
2811 0 : spPSDslp[pp][nn] /= lifetimeTrials;
2812 : }
2813 :
2814 0 : for( size_t nn = 0; nn < sz2Sided; ++nn )
2815 : {
2816 0 : clPSD[nn] =
2817 0 : opdPSD[pp][nn] * norm( ETFlp[pp][nn] ) + pow( sqrtNPSD[pp][nn], 2 ) * norm( NTFlp[pp][nn] );
2818 : }
2819 :
2820 0 : var = sigproc::psdVar( freq2sided, clPSD );
2821 :
2822 0 : pvar = sigproc::psdVar( spFreq, spPSDslp[pp] );
2823 :
2824 0 : for( size_t nn = 0; nn < spPSDslp[pp].size(); ++nn )
2825 : {
2826 0 : spPSDslp[pp][nn] *= var / pvar;
2827 0 : imclp.image( nn )( mnMax + m, mnMax + n ) = spPSDslp[pp][nn];
2828 0 : imclp.image( nn )( mnMax - m, mnMax - n ) = spPSDslp[pp][nn];
2829 : }
2830 : }
2831 :
2832 0 : var = sigproc::psdVar( spFreq, spPSDs[pp] );
2833 :
2834 0 : realT T = ( 1.0 / ( spFreq[1] - spFreq[0] ) ) * 10;
2835 : realT error;
2836 0 : realT tau = pvm( error, spFreq, spPSDs[pp], T ) * ( T ) / var;
2837 0 : taus( mnMax + m, mnMax + n ) = tau;
2838 0 : taus( mnMax - m, mnMax - n ) = tau;
2839 :
2840 0 : var = sigproc::psdVar( spFreq, spPSDslp[pp] );
2841 :
2842 0 : tau = pvm( error, spFreq, spPSDslp[pp], T ) * ( T ) / var;
2843 0 : tauslp( mnMax + m, mnMax + n ) = tau;
2844 0 : tauslp( mnMax - m, mnMax - n ) = tau;
2845 : }
2846 :
2847 : /*********************************************************************/
2848 : // 4.0) Write the results to disk
2849 : /*********************************************************************/
2850 0 : fn = std::format( "{}/speckleLifetimes_{}_si.fits", dir, mags[s] );
2851 : // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2852 0 : ff.write( fn, taus );
2853 :
2854 0 : fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] );
2855 : // dir + "/speckleLifetimes_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2856 0 : ff.write( fn, tauslp );
2857 :
2858 0 : if( writePSDs )
2859 : {
2860 0 : fn = std::format( "{}/specklePSDs_{}_si.fits", dir, mags[s] );
2861 : // dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_si.fits";
2862 0 : ff.write( fn, imc );
2863 :
2864 0 : fn = std::format( "{}/speckleLifetimes_{}_lp.fits", dir, mags[s] );
2865 : // dir + "/specklePSDs_" + ioutils::convert ToString( mags[s] ) + "_lp.fits";
2866 0 : ff.write( fn, imclp );
2867 : }
2868 :
2869 : } // mags
2870 :
2871 0 : return 0;
2872 0 : }
2873 :
2874 : template <typename realT, typename aosysT>
2875 0 : int fourierTemporalPSD<realT, aosysT>::getGridFreq( std::vector<realT> &freq, const std::string &dir )
2876 : {
2877 0 : std::string fn;
2878 0 : fn = dir + "/psds/freq.binv";
2879 0 : return ioutils::readBinVector( freq, fn );
2880 0 : }
2881 :
2882 : template <typename realT, typename aosysT>
2883 0 : int fourierTemporalPSD<realT, aosysT>::getGridPSD( std::vector<realT> &psd, const std::string &dir, int m, int n )
2884 : {
2885 0 : std::string fn;
2886 0 : fn = std::format( "{}/psds/psd_{}_{}.binv", dir, m, n );
2887 : // dir + "/psds/psd_" + ioutils::convert ToString( m ) + '_' + ioutils::convert ToString( n ) + ".binv";
2888 0 : return ioutils::readBinVector( psd, fn );
2889 0 : }
2890 :
2891 : template <typename realT, typename aosysT>
2892 0 : int fourierTemporalPSD<realT, aosysT>::getGridPSD(
2893 : std::vector<realT> &freq, std::vector<realT> &psd, const std::string &dir, int m, int n )
2894 : {
2895 0 : int rv = getGridFreq( freq, dir );
2896 0 : if( rv < 0 )
2897 0 : return rv;
2898 :
2899 0 : return getGridPSD( psd, dir, m, n );
2900 : }
2901 :
2902 : /// Worker function for GSL Integration for the basic sin/cos Fourier modes.
2903 : /** \ingroup mxAOAnalytic
2904 : */
2905 : template <typename realT, typename aosysT>
2906 96450 : realT F_basic( realT kv, void *params )
2907 : {
2908 96450 : fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2909 :
2910 96450 : realT f = Fp->m_f;
2911 96450 : realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2912 :
2913 96450 : realT D = Fp->m_aosys->D();
2914 96450 : realT m = Fp->m_m;
2915 96450 : realT n = Fp->m_n;
2916 96450 : int p = Fp->m_p;
2917 :
2918 96450 : realT ku = f / v_wind;
2919 :
2920 96450 : realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2921 96450 : realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2922 :
2923 96450 : realT Q1 = math::func::jinc( math::pi<realT>() * D * kp );
2924 :
2925 96450 : realT Q2 = math::func::jinc( math::pi<realT>() * D * kpp );
2926 :
2927 96450 : realT Q = ( Q1 + p * Q2 );
2928 :
2929 : realT P =
2930 96450 : Fp->m_aosys->psd( Fp->m_aosys->atm, Fp->_layer_i, sqrt( pow( ku, 2 ) + pow( kv, 2 ) ), Fp->m_aosys->secZeta() );
2931 :
2932 96450 : return P * Q * Q;
2933 : }
2934 :
2935 : template <typename realT>
2936 0 : void turbBoilCubic(
2937 : realT &a, realT &b, realT &c, realT &d, const realT &kv, const realT &f, const realT &Vu, const realT &f0, int pm )
2938 : {
2939 0 : a = Vu * Vu * Vu;
2940 0 : b = -( 3 * Vu * Vu * f + pm * f0 * f0 * f0 );
2941 0 : c = 3 * f * f * Vu;
2942 0 : d = -( f * f * f + pm * f0 * f0 * f0 * kv * kv );
2943 0 : }
2944 :
2945 : /// Worker function for GSL Integration for the modified Fourier modes.
2946 : /** \ingroup mxAOAnalytic
2947 : */
2948 : template <typename realT, typename aosysT>
2949 0 : realT F_mod( realT kv, void *params )
2950 : {
2951 0 : fourierTemporalPSD<realT, aosysT> *Fp = (fourierTemporalPSD<realT, aosysT> *)params;
2952 :
2953 0 : realT f = Fp->m_f;
2954 0 : realT v_wind = Fp->m_aosys->atm.layer_v_wind( Fp->_layer_i );
2955 :
2956 0 : realT D = Fp->m_aosys->D();
2957 0 : realT m = Fp->m_m;
2958 0 : realT n = Fp->m_n;
2959 :
2960 0 : realT f0 = Fp->m_f0;
2961 :
2962 : realT ku;
2963 0 : if( f0 == 0 )
2964 : {
2965 0 : ku = f / v_wind;
2966 :
2967 0 : if( Fp->m_spatialFilter )
2968 : {
2969 : // de-rotate the spatial frequency vector back to pupil coordinates
2970 0 : realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
2971 0 : realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
2972 : // Return if spatially filtered
2973 0 : if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
2974 0 : return 0;
2975 :
2976 0 : if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
2977 0 : return 0;
2978 : }
2979 :
2980 0 : realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
2981 0 : realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
2982 :
2983 0 : realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
2984 :
2985 0 : realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
2986 :
2987 0 : realT QQ = 2 * ( Jp * Jp + Jm * Jm );
2988 :
2989 0 : realT P = Fp->m_aosys->psd( Fp->m_aosys->atm,
2990 0 : Fp->_layer_i,
2991 0 : sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
2992 0 : Fp->m_aosys->lam_sci(),
2993 0 : Fp->m_aosys->lam_wfs(),
2994 0 : Fp->m_aosys->secZeta() );
2995 :
2996 0 : return P * QQ;
2997 : }
2998 : else
2999 : {
3000 : realT a, b, c, d, p, q;
3001 :
3002 0 : turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, 1 );
3003 0 : mx::math::cubicDepressed( p, q, a, b, c, d );
3004 0 : realT t = mx::math::cubicRealRoot( p, q );
3005 :
3006 0 : ku = t - b / ( 3 * a );
3007 :
3008 0 : if( Fp->m_spatialFilter )
3009 : {
3010 : // de-rotate the spatial frequency vector back to pupil coordinates
3011 0 : realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
3012 0 : realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
3013 : // Return if spatially filtered
3014 0 : if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
3015 0 : return 0;
3016 :
3017 0 : if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
3018 0 : return 0;
3019 : }
3020 :
3021 0 : realT kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
3022 0 : realT kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
3023 :
3024 0 : realT Jp = math::func::jinc( math::pi<realT>() * D * kp );
3025 :
3026 0 : realT Jm = math::func::jinc( math::pi<realT>() * D * kpp );
3027 :
3028 0 : realT QQ = 2 * ( Jp * Jp + Jm * Jm );
3029 :
3030 0 : realT P1 = Fp->m_aosys->psd( Fp->m_aosys->atm,
3031 0 : Fp->_layer_i,
3032 0 : sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
3033 0 : Fp->m_aosys->lam_sci(),
3034 0 : Fp->m_aosys->lam_wfs(),
3035 0 : Fp->m_aosys->secZeta() );
3036 :
3037 0 : P1 *= QQ;
3038 :
3039 0 : turbBoilCubic( a, b, c, d, kv, f, v_wind, f0, -1 );
3040 0 : mx::math::cubicDepressed( p, q, a, b, c, d );
3041 0 : t = mx::math::cubicRealRoot( p, q );
3042 :
3043 0 : ku = t - b / ( 3 * a );
3044 :
3045 0 : if( Fp->m_spatialFilter )
3046 : {
3047 : // de-rotate the spatial frequency vector back to pupil coordinates
3048 0 : realT dku = ku * Fp->m_cq - kv * Fp->m_sq;
3049 0 : realT dkv = ku * Fp->m_sq + kv * Fp->m_cq;
3050 : // Return if spatially filtered
3051 0 : if( fabs( dku ) >= Fp->m_aosys->spatialFilter_ku() )
3052 0 : return 0;
3053 :
3054 0 : if( fabs( dkv ) >= Fp->m_aosys->spatialFilter_kv() )
3055 0 : return 0;
3056 : }
3057 :
3058 0 : kp = sqrt( pow( ku + m / D, 2 ) + pow( kv + n / D, 2 ) );
3059 0 : kpp = sqrt( pow( ku - m / D, 2 ) + pow( kv - n / D, 2 ) );
3060 :
3061 0 : Jp = math::func::jinc( math::pi<realT>() * D * kp );
3062 :
3063 0 : Jm = math::func::jinc( math::pi<realT>() * D * kpp );
3064 :
3065 0 : QQ = 2 * ( Jp * Jp + Jm * Jm );
3066 :
3067 0 : realT P2 = Fp->m_aosys->psd( Fp->m_aosys->atm,
3068 0 : Fp->_layer_i,
3069 0 : sqrt( pow( ku, 2 ) + pow( kv, 2 ) ),
3070 0 : Fp->m_aosys->lam_sci(),
3071 0 : Fp->m_aosys->lam_wfs(),
3072 0 : Fp->m_aosys->secZeta() );
3073 :
3074 0 : P2 *= QQ;
3075 :
3076 0 : return 0.5 * ( P1 + P2 );
3077 : }
3078 : }
3079 :
3080 : /*extern template
3081 : struct fourierTemporalPSD<float, aoSystem<float, vonKarmanSpectrum<float>, std::ostream>>;*/
3082 :
3083 : extern template struct fourierTemporalPSD<double, aoSystem<double, vonKarmanSpectrum<double>, std::ostream>>;
3084 :
3085 : /*
3086 : extern template
3087 : struct fourierTemporalPSD<long double, aoSystem<long double, vonKarmanSpectrum<long double>, std::ostream>>;
3088 :
3089 : #ifdef HASQUAD
3090 : extern template
3091 : struct fourierTemporalPSD<__float128, aoSystem<__float128, vonKarmanSpectrum<__float128>, std::ostream>>;
3092 : #endif
3093 : */
3094 :
3095 : } // namespace analysis
3096 : } // namespace AO
3097 : } // namespace mx
3098 :
3099 : #endif // fourierTemporalPSD_hpp
|
