mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
simulatedAOSystem.hpp
Go to the documentation of this file.
1 /** \file
2  * \author Jared R. Males (jaredmales@gmail.com)
3  * \brief
4  * \ingroup mxAO_sim_files
5  *
6  */
7 
8 #ifndef simulatedAOSystem_hpp
9 #define simulatedAOSystem_hpp
10 
11 #include <iostream>
12 #include <fstream>
13 
14 #include <cuda_runtime.h>
15 #include <cublas_v2.h>
16 
17 
18 #include "../../improc/imagePads.hpp"
19 #include "../../wfp/imagingUtils.hpp"
20 #include "../../wfp/fraunhoferPropagator.hpp"
21 
22 #include "../../ioutils/fits/fitsFile.hpp"
23 #include "../../ioutils/fits/fitsUtils.hpp"
24 
25 #include "../../improc/eigenImage.hpp"
26 #include "../../improc/eigenCube.hpp"
27 
28 #include "../../sys/timeUtils.hpp"
29 #include "../../sigproc/signalWindows.hpp"
30 
31 #include "wavefront.hpp"
32 #include "../aoPaths.hpp"
33 
34 
35 #ifdef DEBUG
36 #define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
37 #else
38 #define BREAD_CRUMB
39 #endif
40 
41 namespace mx
42 {
43 namespace AO
44 {
45 namespace sim
46 {
47 
48 
49 
50 ///A simulated AO system.
51 /**
52  *
53  * Minimum requirement for _turbSeqT:
54  \code
55  {
56  int wfPS(realT ps); //Sets the wavefront platescale (only needed for cascaded systems).
57  int frames(); //Returns the number of frames, sets the number of iterations to run.
58  int nextWF(wavefront<realT> & wf); //Fills in the wavefront with phase and amplitude
59  };
60  \endcode
61  *
62  */
63 template<typename _realT, typename _wfsT, typename _reconT, typename _filterT, typename _dmT, typename _turbSeqT, typename _coronT>
64 class simulatedAOSystem
65 {
66 public:
67 
68  typedef _realT realT; ///< The floating point type used for all calculations
69 
70  typedef mx::AO::sim::wavefront<realT> wavefrontT; ///< The wavefront type
71 
72  typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT; ///<The real image type
73 
74  typedef wfp::imagingArray<std::complex<realT>, wfp::fftwAllocator<std::complex<realT> >, 0> complexImageT; ///<The complex image type
75 
76  typedef _wfsT wfsT;
77  typedef _reconT reconT;
78  typedef _filterT filterT;
79  typedef _dmT dmT;
80  typedef _turbSeqT turbSeqT;
81  typedef _coronT coronT;
82 
83  //cuda admin
84  cublasHandle_t m_cublasHandle;
85 
86 
87  std::string _sysName; ///< The system name for use in mx::AO::path
88  std::string _wfsName; ///< The WFS name for use in the mx::AO::path
89  std::string _pupilName; ///< The pupil name for use in the mx::AO::path
90 
91  bool m_sfImagePlane {false};
92 
93  wfsT wfs;
94  reconT recon;
95  filterT filter;
96  dmT dm;
97  turbSeqT turbSeq;
98 
99 
100 
101 
102 
103  void wfPS(realT wps)
104  {
105  }
106 
107 
108 
109  realT _wfsLambda;
110  realT _sciLambda;
111 
112  long _frameCounter;
113 
114  bool _loopClosed;
115  int _loopClosedDelay;
116  int _lowOrdersDelay;
117 
118  std::string _rmsFile;
119  std::ofstream _rmsOut;
120  bool _rmsUnwrap;
121 
122  std::string _ampFile;
123  //std::ofstream m_ampOut;
124 
125  //Members which have implemented accessors are moved here as I go:
126 protected:
127  realT m_simStep; ///< The simulation step size in seconds.
128 
129  realT m_D;
130 
131  realT m_wfPS;
132 
133  int m_wfSz;
134 
135 
136  mx::improc::eigenImage<realT> m_ampOut;
137 
138 public:
139  /** Wavefront Outputs
140  * @{
141  */
142 
143  bool m_writeWavefronts {true};
144  std::string m_wfFileBase {"simAOWF"};
145  mx::improc::eigenCube<realT> m_wfPhase;
146  mx::improc::eigenCube<realT> m_wfAmp;
147 
148  int m_nWFPerFile {500};
149  int m_currWF {0};
150  int m_currWFFile {0};
151 
152  ///@}
153 
154  /** Image outputs
155  * @{
156  */
157 
158  int m_saveSz {0};
159  bool _writeIndFrames;
160 
161  long _saveFrameStart;
162 
163  std::string _psfFileBase;
164  mx::improc::eigenCube<realT> _psfs;
165  imageT _psfOut;
166 
167 
168  mx::improc::eigenCube<realT> _corons;
169  imageT _coronOut;
170 
171  int _nPerFile;
172  int _currImage;
173 
174  int _currFile;
175 
176  complexImageT _complexPupil;
177  complexImageT _complexPupilCoron;
178  //complexImageT _complexFocal;
179  imageT _realFocal;
180 
181  imageT _realPupil;
182  imageT _realFocalCoron;
183 
184  coronT m_coronagraph;
185 
186  //wfp::fraunhoferPropagator<complexImageT> _fi;
187 
188 
189  std::vector<typename filterT::commandT> _delayedCommands;
190  int _commandDelay;
191  std::vector<int> _goodCommands;
192 
193  /** @}
194  */
195 
196  ///Default c'tor
197  simulatedAOSystem();
198 
199  ///Destructor
200  ~simulatedAOSystem();
201 
202 
203  ///Initialize the system
204  /**
205  * Initializes the basic parts of the system.
206  *
207  * \returns 0 on success
208  * \returns -1 on an error (simulation should not continue if this happens).
209  *
210  */
211  int initSystem( const std::string & sysName, ///< [in] System name.
212  typename dmT::specT & dmSpec, ///< [in] DM Specification.
213  const std::string & wfsName, ///< [in] WFS Name.
214  const std::string & pupilName, ///< [in] Name of the system pupil.
215  const int & wfSz ///< [in] Size of the wavefront used for propagation.
216  );
217 
218  void initSim( typename reconT::specT & reconSpec,
219  realT simStep,
220  int commandDelay,
221  const std::string & coronName
222  );
223 
224  imageT _pupil; ///< The system pupil. This is generally a binary mask and will be applied at various points in the propagation.
225  imageT _pupilMask; ///< A pupil mask which is applied once at the beginning of propagation. Could be apodized and/or different from _pupil.
226 
227  imageT _postMask;
228  imageT _coronPhase;
229 
230 
231  int _npix;
232 
233  ///Measure the system response matrix
234  /** System should be initialized with initSystemCal.
235  *
236  * \param amp
237  * \param rmatName
238  * \param nmodes
239  */
240  void takeResponseMatrix(realT amp, std::string rmatName, int nmodes=0);
241 
242  int frames();
243 
244  void calcOpenLoopAmps(wavefrontT & wf);
245 
246  void nextWF(wavefrontT & wf);
247 
248  void runTurbulence();
249 
250  /** \name Member Access
251  * @{
252  */
253 
254  /// Set the simulation step size.
255  /** Units of step size are seconds.
256  *
257  * \returns 0 on succces
258  * \returns -1 on error [currently none]
259  */
260  int simStep( const double & ss /**< [in] the new value of simulation stepsize [sec] */);
261 
262  /// Get the simulation step size.
263  /**
264  * \returns the current value of m_simStep [sec].
265  */
266  double simStep();
267 
268  int D( const realT & nD );
269 
270  realT D();
271 
272  int wfPS (const realT & nwfPS);
273 
274  realT wfPS();
275 
276  int Dpix();
277 
278  int wfSz( const int & nwfSz);
279 
280  int wfSz();
281 
282  bool m_doCoron {false};
283 
284  ///@}
285 
286 public:
287  double dt_turbulence {0};
288  double dt_wfs {0};
289  double dt_recon {0};
290  double dt_dmcomb {0};
291  double dt_total {0};
292 
293 };
294 
295 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
296 simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::simulatedAOSystem()
297 {
298  _frameCounter = 0;
299  _loopClosed = false;
300  _loopClosedDelay = 0;
301  _lowOrdersDelay = 0;
302 
303  _rmsUnwrap = false;
304 
305  _writeIndFrames = false;
306  _saveFrameStart = 0;
307  _nPerFile = 100;
308  _currImage = 0;
309  _currFile = 0;
310 
311  m_coronagraph.m_fileDir = sys::getEnv("MX_AO_DATADIR") + "/" + "coron/";
312 
313  _npix = 0;
314 
315 
316  cublasCreate(&m_cublasHandle);
317  cublasSetPointerMode(m_cublasHandle, CUBLAS_POINTER_MODE_DEVICE);//CUBLAS_POINTER_MODE_HOST
318 }
319 
320 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
321 simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::~simulatedAOSystem()
322 {
323  if(m_ampOut.rows() > 0 && _ampFile != "")
324  {
325  mx::fits::fitsFile<realT> ff;
326  ff.write(_ampFile+".fits", m_ampOut);
327  }
328 
329  if(_rmsOut.is_open()) _rmsOut.close();
330 #if 1
331  //Write out the psf and coron cubes one last time
332  if(_psfFileBase !="" && _currImage > 0 && _writeIndFrames)
333  {
334  mx::fits::fitsFile<realT> ff;
335 
336  std::string fn = _psfFileBase + "_psf_" + ioutils::convertToString(_currFile) + ".fits";
337 
338  BREAD_CRUMB;
339 
340  ff.write(fn, _psfs);//_psfs.data(), _psfs.rows(), _psfs.cols(), _currImage);
341 
342  if(m_doCoron)
343  {
344  std::string fn = _psfFileBase + "_coron_" + ioutils::convertToString(_currFile) + ".fits";
345 
346  BREAD_CRUMB;
347 
348  ff.write(fn, _corons);// _corons.data(), _corons.rows(), _corons.cols(), _currImage);
349  }
350  }
351 #endif
352 
353 }
354 
355 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
356 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::initSystem( const std::string & sysName,
357  typename dmT::specT & dmSpec,
358  const std::string & wfsName,
359  const std::string & pupilName,
360  const int & wfSz )
361 {
362  _sysName = sysName;
363  _wfsName = wfsName;
364  _pupilName = pupilName;
365 
366  fits::fitsFile<realT> ff;
367  fits::fitsHeader head;
368 
369  //Initialize the pupil.
370 
371  std::string pupilFile = mx::AO::path::pupil::pupilFile(_pupilName);
372 
373  ff.read(_pupil, head, pupilFile);
374 
375  m_D = head["PUPILD"].Value<realT>(); //pupilD;
376  m_wfPS = head["SCALE"].Value<realT>();
377 
378  //Set the wavefront size.
379  m_wfSz = wfSz;
380  wfs.wfSz(m_wfSz);
381 
382  //Turbulence sequence
383  turbSeq._pupil = &_pupil;
384  turbSeq.wfPS(m_wfPS);
385 
386  //DM Initialization.
387  dm.initialize( *this, dmSpec, _pupilName);
388 
389  wfs.linkSystem(*this);
390 
391  _loopClosed = false;
392 
393  return 0;
394 
395 }
396 
397 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
398 void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::initSim( typename reconT::specT & reconSpec,
399  realT simStep,
400  int commandDelay,
401  const std::string & coronName
402  )
403 {
404 
405  fits::fitsFile<realT> ff;
406  fits::fitsHeader head;
407 
408  m_simStep = simStep;
409  wfs.simStep(m_simStep);
410 
411 
412  filter.initialize(dm.nModes());
413 
414  recon.initialize(*this, reconSpec);
415  dm.calAmp(recon.calAmp());
416 
417  _loopClosed = false;
418 
419  _commandDelay = commandDelay;
420 
421  _delayedCommands.resize((_commandDelay+1)*5);
422  _goodCommands.resize((_commandDelay+1)*5);
423 
424  if(coronName != "")
425  {
426  m_coronagraph.wfSz(m_wfSz);
427  m_coronagraph.loadCoronagraph(coronName);
428  }
429 }
430 
431 
432 // template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
433 // void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::initSystemCal( const std::string & sysName,
434 // const std::string & dmName,
435 // const std::string & wfsName,
436 // const std::string & pupilName,
437 // const std::string & basisName,
438 // const bool & basisOrtho,
439 // const int & wfSz )
440 // {
441 // _sysName = sysName;
442 // _dmName = dmName;
443 // _wfsName = wfsName;
444 // _pupilName = pupilName;
445 // _basisName = basisName;
446 // _basisOrtho = basisOrtho;
447 //
448 // mx::fitsFile<realT> ff;
449 // mx::fitsHeader head;
450 //
451 //
452 // std::string pupilFile = mx::AO::path::pupil::pupilFile(_pupilName);
453 //
454 // ff.read(pupilFile, _pupil, head);
455 //
456 // //m_D = head["SCALE"].Value<realT>(); //pupilD;
457 // m_D = head["PUPILD"].Value<realT>(); //pupilD;
458 // m_wfPS = head["SCALE"].Value<realT>();
459 //
460 // m_wfSz = wfSz;
461 // wfs.wfSz(m_wfSz);
462 //
463 // m_wfPS = m_D/std::max(_pupil.rows(), _pupil.cols());
464 // turbSeq.wfPS(m_wfPS);
465 //
466 // //DM Initialization.
467 // dm.initialize( _dmName, _basisName, _basisOrtho, _pupilName);
468 //
469 // //dm.loadModes(basisSet, pupil);
470 //
471 // wfs.linkSystem(*this);
472 //
473 // _loopClosed = false;
474 //
475 //
476 // }
477 
478 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
479 void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::takeResponseMatrix( realT amp,
480  std::string rmatID,
481  int nmodes )
482 {
483  complexImageT cpup;
484  cpup.resize(m_wfSz, m_wfSz);
485 
486  recon.initializeRMat(dm.nModes(), amp, wfs.detRows(), wfs.detCols());
487 
488  typename filterT::commandT measureVec;
489 
490  wavefrontT currWF;
491  currWF.setAmplitude(_pupil);
492 
493  wfs.iTime(1);
494  wfs.roTime(1);
495  wfs.detector.noNoise(true);
496 
497  double tO,tF, t0, t1, t_applyMode = 0, t_senseWF = 0, t_calcMeas = 0, t_accum = 0;
498 
499  std::cerr << dm.nModes() << "\n";
500 
501  tO = sys::get_curr_time();
502 
503  for(int i=0;i< dm.nModes(); ++i)
504  {
505  BREAD_CRUMB;
506 
507  currWF.setPhase(_pupil*0);
508 
509  realT s_amp = amp;
510  if(nmodes>0 && i>= nmodes)
511  {
512  s_amp =0;
513  wfs.detectorImage.image.setZero();
514  }
515  else
516  {
517  t0 = sys::get_curr_time();
518  dm.applyMode(currWF, i, s_amp, 0.8e-6);
519  t1 = sys::get_curr_time();
520  t_applyMode += t1-t0;
521 
522  BREAD_CRUMB;
523 
524  t0 = sys::get_curr_time();
525  wfs.senseWavefrontCal(currWF);
526  t1 = sys::get_curr_time();
527 
528  t_senseWF += t1-t0;
529  }
530 
531  BREAD_CRUMB;
532 
533  t0 = sys::get_curr_time();
534 
535  recon.calcMeasurement(measureVec, wfs.detectorImage);
536 
537  t1 = sys::get_curr_time();
538 
539  t_calcMeas += t1-t0;
540 
541  BREAD_CRUMB;
542 
543  t0 = sys::get_curr_time();
544  recon.accumulateRMat(i, measureVec, wfs.detectorImage);
545  t1 = sys::get_curr_time();
546  t_accum += t1-t0;
547 
548  BREAD_CRUMB;
549 
550  }
551  tF = sys::get_curr_time();
552  std::cout << ( (realT) dm.nModes())/(tF - tO) << " Hz\n";
553 
554  std::cout << t_applyMode/dm.nModes() << " " << t_senseWF/dm.nModes() << " " << t_calcMeas/dm.nModes() << " " << t_accum/dm.nModes();
555  std::cout << " " << dm.t_mm/dm.nModes() << " " << dm.t_sum/dm.nModes() << " " << "\n";
556 
557  std::string fname;
558  fname = mx::AO::path::sys::cal::rMat(_sysName, dm.name(), _wfsName, _pupilName, dm.basisName(), rmatID, true);
559 
560  std::cout << fname << "\n";
561  recon.saveRMat(fname);
562 
563  fname = mx::AO::path::sys::cal::rImages(_sysName, dm.name(), _wfsName, _pupilName, dm.basisName(), rmatID, true);
564  recon.saveRImages(fname);
565 
566 }
567 
568 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
569 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::frames()
570 {
571  return turbSeq.frames();
572 }
573 
574 /*
575 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
576 void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::calcOpenLoopAmps(wavefrontT & wf)
577 {
578  BREAD_CRUMB;
579 
580  int npix = _pupil->sum();
581 
582  imageT olAmps.resize(1, dm.nModes());
583 
584  #pragma omp parallel for
585  for(int j=0; j< dm.nModes; ++j)
586  {
587  realT amp;
588 
589 
590  amp = (wf.phase*dm._infF->image(j)).sum()/ npix;
591 
592  olAmps(0,j) = amp;
593  }
594 
595 
596 }
597 */
598 
599 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
600 void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::nextWF(wavefrontT & wf)
601 {
602 
603  realT rms_ol, rms_cl;
604 
605  BREAD_CRUMB;
606 
607  if(_npix == 0)
608  {
609  _npix = _pupil.sum();
610  }
611 
612  BREAD_CRUMB;
613 
614  double t0 = sys::get_curr_time();
615  turbSeq.nextWF(wf);
616  dt_turbulence += sys::get_curr_time() - t0;
617 
618  wf.iterNo = _frameCounter;
619 
620  BREAD_CRUMB;
621 
622  //Mean subtraction on the system pupil.
623  realT mn = (wf.phase * _pupil).sum()/_npix;
624 
625  //Apply the pupil mask just once.
626  wf.phase = (wf.phase-mn)*_pupil;
627 
628  if(_pupilMask.rows() == _pupil.rows() && _pupilMask.cols() == _pupil.cols())
629  {
630  wf.phase *= _pupilMask;
631  }
632 
633 
634  BREAD_CRUMB;
635 
636  if(_frameCounter % 10 == 0) rms_ol = sqrt( wf.phase.square().sum()/ _npix );
637 
638  BREAD_CRUMB;
639 
640  typename filterT::commandT measuredAmps, commandAmps;
641 
642  filter.initMeasurements(measuredAmps, commandAmps);
643 
644  if(_frameCounter == 0)
645  {
646  for(size_t i=0; i<_delayedCommands.size();++i)
647  {
648  _goodCommands[i] = 0;
649  }
650  }
651 
652  BREAD_CRUMB;
653  if(_loopClosed)
654  {
655  dm.applyShape(wf, _wfsLambda);
656 
657  bool newCV;
658 
659  BREAD_CRUMB;
660 
661  t0 = sys::get_curr_time();
662  newCV = wfs.senseWavefront(wf);
663  dt_wfs += sys::get_curr_time() - t0;
664 
665  if(newCV)
666  {
667  BREAD_CRUMB;
668 
669  t0 = sys::get_curr_time();
670  recon.reconstruct(measuredAmps, wfs.m_detectorImage);
671  dt_recon += sys::get_curr_time() - t0;
672 
673  BREAD_CRUMB;
674 
675  //Record amps if we're saving
676  if(_ampFile != "")
677  {
678  //Check if initialized
679  if(m_ampOut.cols() != turbSeq.frames() || m_ampOut.rows() != measuredAmps.measurement.size()+1)
680  {
681  if(m_ampOut.cols() != 0)
682  {
683  std::cerr << "m_ampOut already allocated but cols not frames\n";
684  exit(0);
685  }
686 
687  if(m_ampOut.rows() != 0)
688  {
689  std::cerr << "m_ampOut already allocated but rows not measurements\n";
690  exit(0);
691  }
692 
693  m_ampOut.resize(measuredAmps.measurement.size()+1, turbSeq.frames());
694  m_ampOut.setZero();
695  }
696 
697  int n = floor(measuredAmps.iterNo);
698  m_ampOut(0,n) = measuredAmps.iterNo;
699  for(int i=1;i<measuredAmps.measurement.size()+1; ++i)
700  {
701  m_ampOut(i,n) = measuredAmps.measurement[i-1];
702  }
703  }
704 
705  BREAD_CRUMB;
706 
707  int nAmps = ( _frameCounter % _delayedCommands.size() );
708 
709  _delayedCommands[nAmps].measurement = measuredAmps.measurement;
710  _delayedCommands[nAmps].iterNo = measuredAmps.iterNo;
711 
712  _goodCommands[nAmps] = 1;
713  }
714  else
715  {
716  int nAmps = ( _frameCounter % _delayedCommands.size() );
717  _goodCommands[nAmps] = 0;
718  }
719 
720 
721  int _currCommand = ( _frameCounter % _delayedCommands.size() ) - _commandDelay;
722  if(_currCommand < 0) _currCommand += _delayedCommands.size();
723 
724  if(_goodCommands[_currCommand])
725  {
726  BREAD_CRUMB;
727 
728  filter.filterCommands(commandAmps, _delayedCommands[_currCommand], _frameCounter);
729 
730  BREAD_CRUMB;
731 
732  t0=sys::get_curr_time();
733  dm.setShape(commandAmps);
734  dt_dmcomb += sys::get_curr_time() - t0;
735  }
736 
737 
738  }
739  else
740  {
741  int nAmps = ( _frameCounter % _delayedCommands.size() );
742  _goodCommands[nAmps] = 0;
743  }
744 
745 
746  //**** If the _postMask isn't set, set it to _pupil ****//
747  if(_postMask.rows() != _pupil.rows())
748  {
749  _postMask = _pupil;
750  }
751 
752 
753  //**** Calculate RMS phase ****//
754  if(_frameCounter % 10 == 0)
755  {
756  mn = (wf.phase * _pupil).sum()/_npix;
757  rms_cl = sqrt( (wf.phase-mn).square().sum()/ _postMask.sum() );
758  std::cout << _frameCounter << " WFE: " << rms_ol << " " << rms_cl << " [rad rms phase]\n";
759  }
760 
761  if(m_sfImagePlane)
762  {
763  wfs.filter().filter(wf.phase);
764  }
765 
766  BREAD_CRUMB;
767 
768  if(_rmsFile != "")
769  {
770  if(! _rmsOut.is_open() )
771  {
772  _rmsOut.open(_rmsFile);
773  _rmsOut << "#open-loop-wfe closed-loop-wfe [rad rms phase]\n";
774  }
775  _rmsOut << rms_ol << " " << rms_cl << std::endl;
776  }
777 
778  BREAD_CRUMB;
779 
780  if( m_writeWavefronts )
781  {
782  if( m_wfPhase.rows() != wf.phase.rows() || m_wfPhase.cols() != wf.phase.cols() || m_wfPhase.planes() != m_nWFPerFile ||
783  m_wfAmp.rows() != wf.amplitude.rows() || m_wfAmp.cols() != wf.amplitude.cols() || m_wfAmp.planes() != m_nWFPerFile )
784  {
785  m_wfPhase.resize( wf.phase.rows(), wf.phase.cols(), m_nWFPerFile);
786  m_wfAmp.resize( wf.amplitude.rows(), wf.amplitude.cols(), m_nWFPerFile);
787 
788  m_currWF = 0;
789  }
790 
791  m_wfPhase.image(m_currWF) = wf.phase;
792  m_wfAmp.image(m_currWF) = wf.amplitude;
793 
794  ++m_currWF;
795 
796  if( m_currWF >= m_nWFPerFile )
797  {
798  std::cerr << "Write to WF file here . . . \n";
799 
800  mx::fits::fitsFile<realT> ff;
801 
802  std::string fn = m_wfFileBase + "_phase_" + ioutils::convertToString<int,5,'0'>(m_currWFFile) + ".fits";
803  ff.write(fn, m_wfPhase);
804 
805  fn = m_wfFileBase + "_amp_" + ioutils::convertToString<int,5,'0'>(m_currWFFile) + ".fits";
806  //ff.write(fn, m_wfAmp);
807 
808  ++m_currWFFile;
809  m_currWF = 0;
810  }
811 
812  }
813 
814  BREAD_CRUMB;
815 
816  if(_psfFileBase != "" && _frameCounter > _saveFrameStart)
817  {
818  if(_psfOut.rows() == 0 || (_psfs.planes() == 0 && _writeIndFrames))
819  {
820  if(m_saveSz <= 0) m_saveSz = m_wfSz;
821 
822  if(_writeIndFrames) _psfs.resize(m_saveSz, m_saveSz, _nPerFile);
823  _psfOut.resize(m_saveSz, m_saveSz);
824  _psfOut.setZero();
825 
826  _realFocal.resize(m_saveSz, m_saveSz);
827 
828  if(m_doCoron)
829  {
830  if(_writeIndFrames) _corons.resize(m_saveSz, m_saveSz, _nPerFile);
831  _coronOut.resize(m_saveSz, m_saveSz);
832  _coronOut.setZero();
833 
834  _realFocalCoron.resize(m_saveSz, m_saveSz);
835  }
836  }
837 
838  BREAD_CRUMB;
839 
840  //Propagate Coronagraph
841  if(m_doCoron)
842  {
843  wf.getWavefront(_complexPupilCoron, m_wfSz);
844 
845  m_coronagraph.propagate(_realFocalCoron, _complexPupilCoron);
846 
847  if(_writeIndFrames) _corons.image(_currImage) = _realFocalCoron;
848 
849  _coronOut += _realFocalCoron;
850  }
851 
852  //Propagate PSF
853  wf.getWavefront(_complexPupil, m_wfSz);
854 
855  m_coronagraph.propagateNC(_realFocal, _complexPupil);
856 
857  if(_writeIndFrames) _psfs.image(_currImage) = _realFocal;
858 
859  _psfOut += _realFocal;
860 
861 
862 
863  BREAD_CRUMB;
864 
865 
866 
867 
868  ++_currImage;
869 #if 1
870  if(_currImage >= _nPerFile && _writeIndFrames)
871  {
872  mx::fits::fitsFile<realT> ff;
873 
874  std::string fn = _psfFileBase + "_psf_" + ioutils::convertToString(_currFile) + ".fits";
875 
876  BREAD_CRUMB;
877 
878  ff.write(fn, _psfs);// _psfs.data(), _psfs.rows(), _psfs.cols(), _psfs.planes());
879 
880  if(m_doCoron)
881  {
882  std::string fn = _psfFileBase + "_coron_" + ioutils::convertToString(_currFile) + ".fits";
883 
884  BREAD_CRUMB;
885 
886  ff.write(fn, _corons);//_corons.data(), _corons.rows(), _corons.cols(), _corons.planes());
887  }
888 
889  ++_currFile;
890  _currImage = 0;
891  }
892 #endif
893  }//if(_psfFileBase != "" && _frameCounter > _saveFrameStart)
894 
895  ++_frameCounter;
896 
897 }//void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::nextWF(wavefrontT & wf)
898 
899 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
900 void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::runTurbulence()
901 {
902  wavefrontT currWF;
903 
904  if(m_doCoron)
905  {
906  if(m_coronagraph.wfSz() != m_wfSz)
907  {
908  std::cerr << "doCoron is true, but coronagraph not loaded or wavefront sizes don't match.\n";
909  std::cerr << "Continuing . . .\n";
910  }
911  }
912 
913 
914  _wfsLambda = wfs.lambda();
915 
916  double t0 = sys::get_curr_time();
917  for(size_t i=0;i<turbSeq.frames();++i)
918  {
919  //std::cout << i << "/" << turbSeq.frames() << "\n" ;
920 
921 
922  if(i == 0)
923  {
924  turbSeq._loopClosed = true;
925  //wfs.applyFilter = false;
926  }
927 
928  if(i == (size_t) _loopClosedDelay) _loopClosed = true;
929 
930 
931 
932  BREAD_CRUMB;
933 
934  nextWF(currWF);
935 
936  BREAD_CRUMB;
937 
938  }
939 
940  double dt_total = sys::get_curr_time() - t0;
941 
942  if(_psfFileBase != "")
943  {
944  fits::fitsFile<realT> ff;
945  std::string fn = _psfFileBase + "_psf.fits";
946  BREAD_CRUMB;
947  ff.write(fn, _psfOut);
948 
949  if(m_doCoron)
950  {
951  BREAD_CRUMB;
952  fn = _psfFileBase + "_coron.fits";
953  ff.write(fn, _coronOut);
954  }
955  }
956 
957  BREAD_CRUMB;
958 
959  std::cout << "Timing: \n";
960  std::cout << " Total Time: " << dt_total << " sec / " << turbSeq.frames()/dt_total << " fps\n";
961  std::cout << " Turbulence: " << dt_turbulence << " sec / " << turbSeq.frames()/dt_turbulence << " fps / " << dt_turbulence/dt_total << "%\n";
962  std::cout << " WFS: " << dt_wfs << " sec / " << turbSeq.frames()/dt_wfs << " fps / " << dt_wfs/dt_total << "%\n";
963  std::cout << " Recon: " << dt_recon << " sec / " << turbSeq.frames()/dt_recon << " fps / " << dt_recon/dt_total << "%\n";
964  std::cout << " dmcomb: " << dt_dmcomb << " sec / " << turbSeq.frames()/dt_dmcomb << " fps / " << dt_dmcomb/dt_total << "%\n";
965 
966 
967 }//void simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::runTurbulence()
968 
969 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
970 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::simStep( const double & ss)
971 {
972  m_simStep = ss;
973 
974  return 0;
975 }
976 
977 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
978 double simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::simStep()
979 {
980  return m_simStep;
981 }
982 
983 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
984 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::D( const realT & nD )
985 {
986  m_D = nD;
987 
988  return 0;
989 }
990 
991 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
992 realT simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::D()
993 {
994  return m_D;
995 }
996 
997 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
998 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::wfPS (const realT & nwfPS)
999 {
1000  m_wfPS = nwfPS;
1001 
1002  return 0;
1003 }
1004 
1005 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
1006 realT simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::wfPS()
1007 {
1008  return m_wfPS;
1009 }
1010 
1011 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
1012 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::Dpix()
1013 {
1014  return (int)(m_D/m_wfPS +std::numeric_limits<realT>::epsilon());
1015 }
1016 
1017 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
1018 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::wfSz( const int & nwfSz)
1019 {
1020  m_wfSz = nwfSz;
1021 
1022  return 0;
1023 }
1024 
1025 template<typename realT, typename wfsT, typename reconT, typename filterT, typename dmT, typename turbSeqT, typename coronT>
1026 int simulatedAOSystem<realT, wfsT, reconT, filterT, dmT, turbSeqT, coronT>::wfSz()
1027 {
1028  return m_wfSz;
1029 }
1030 
1031 } //namespace sim
1032 } //namespace AO
1033 } //namespace mx
1034 
1035 #endif //simulatedAOSystem_hpp
mx::AO::sim::simulatedAOSystem::initSystem
int initSystem(const std::string &sysName, typename dmT::specT &dmSpec, const std::string &wfsName, const std::string &pupilName, const int &wfSz)
Initialize the system.
Definition: simulatedAOSystem.hpp:356
mx::AO::sim::simulatedAOSystem::_wfsName
std::string _wfsName
The WFS name for use in the mx::AO::path.
Definition: simulatedAOSystem.hpp:88
mx::AO::sim::simulatedAOSystem::_pupilName
std::string _pupilName
The pupil name for use in the mx::AO::path.
Definition: simulatedAOSystem.hpp:89
mx::AO::sim::simulatedAOSystem::realT
_realT realT
The floating point type used for all calculations.
Definition: simulatedAOSystem.hpp:68
mx::AO::sim::simulatedAOSystem::m_simStep
realT m_simStep
The simulation step size in seconds.
Definition: simulatedAOSystem.hpp:127
mx::AO::sim::simulatedAOSystem::imageT
Eigen::Array< realT, Eigen::Dynamic, Eigen::Dynamic > imageT
The real image type.
Definition: simulatedAOSystem.hpp:72
mx::AO::sim::simulatedAOSystem::_pupilMask
imageT _pupilMask
A pupil mask which is applied once at the beginning of propagation. Could be apodized and/or differen...
Definition: simulatedAOSystem.hpp:225
mx::AO::sim::simulatedAOSystem::_sysName
std::string _sysName
The system name for use in mx::AO::path.
Definition: simulatedAOSystem.hpp:87
mx::AO::sim::simulatedAOSystem::_pupil
imageT _pupil
The system pupil. This is generally a binary mask and will be applied at various points in the propag...
Definition: simulatedAOSystem.hpp:224
mx::AO::sim::simulatedAOSystem::simStep
double simStep()
Get the simulation step size.
Definition: simulatedAOSystem.hpp:978
mx::AO::sim::simulatedAOSystem::takeResponseMatrix
void takeResponseMatrix(realT amp, std::string rmatName, int nmodes=0)
Measure the system response matrix.
Definition: simulatedAOSystem.hpp:479
mx::AO::sim::simulatedAOSystem::complexImageT
wfp::imagingArray< std::complex< realT >, wfp::fftwAllocator< std::complex< realT > >, 0 > complexImageT
The complex image type.
Definition: simulatedAOSystem.hpp:74
mx::AO::sim::simulatedAOSystem::wavefrontT
mx::AO::sim::wavefront< realT > wavefrontT
The wavefront type.
Definition: simulatedAOSystem.hpp:70
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::fits::fitsHeader
Class to manage a complete fits header.
Definition: fitsHeader.hpp:52
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::convertToString
std::string convertToString(const typeT &value, int precision=0)
Convert a numerical value to a string.
Definition: stringUtils.hpp:82
mx::sys::getEnv
std::string getEnv(const std::string &estr)
Return the value of an environment variable.
Definition: environment.cpp:33
mx::sys::get_curr_time
typeT get_curr_time(timespec &tsp)
Get the current system time in seconds.
Definition: timeUtils.hpp:63
mx
The mxlib c++ namespace.
Definition: mxError.hpp:107
mx::AO::sim::wavefront
Structure containing the phase and amplitude of a wavefront
Definition: wavefront.hpp:26