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