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