Line data Source code
1 : /** \file turbSequence.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Manage a turbulence sequence saved as FITS files.
4 : * \ingroup mxAO_sim_files
5 : *
6 : */
7 :
8 : #ifndef turbSequence_hpp
9 : #define turbSequence_hpp
10 :
11 : #include <vector>
12 : #include <string>
13 : #include <iostream>
14 :
15 : #include <Eigen/Dense>
16 :
17 : #include "../../ioutils/fileUtils.hpp"
18 : #include "../../ioutils/fits/fitsFile.hpp"
19 : #include "../../error/mxErrorOld.hpp"
20 : #include "../../improc/eigenCube.hpp"
21 :
22 : #include "wavefront.hpp"
23 :
24 : namespace mx
25 : {
26 : namespace AO
27 : {
28 : namespace sim
29 : {
30 :
31 : template <typename _realT>
32 : struct turbSequence
33 : {
34 : typedef _realT realT;
35 : typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
36 :
37 : std::vector<std::string> _phaseFnames;
38 : std::vector<std::string> _ampFnames;
39 :
40 : int _size{ 0 };
41 :
42 : int _nPerCube{ 0 };
43 :
44 : int _files{ 0 };
45 :
46 : int _frames{ 0 };
47 :
48 : int _currFileNo{ 0 };
49 : int _currFrameNo{ 0 };
50 :
51 : realT _wfPS{ 1 };
52 : realT _F0Photons{ 2e11 };
53 :
54 : realT _starMag{ 0 };
55 :
56 : realT _pixVal;
57 :
58 : imageT *_pupil;
59 :
60 : bool _loopClosed;
61 :
62 : improc::eigenCube<realT> _currPhase;
63 : improc::eigenCube<realT> _currAmp;
64 :
65 : bool _phaseOnly{ false };
66 :
67 0 : turbSequence()
68 0 : {
69 0 : }
70 :
71 : /// Return the size of the turbulence screens
72 : /**
73 : * \returns _size
74 : */
75 0 : int size()
76 : {
77 0 : return _size;
78 : }
79 :
80 : /// Return the number of frames per cube
81 : /**
82 : * \returns _nPerCube
83 : */
84 0 : int nPerCube()
85 : {
86 0 : return _nPerCube;
87 : }
88 :
89 : /// Return the number of files
90 : /**
91 : * \returns _files
92 : */
93 0 : int files()
94 : {
95 0 : return _files;
96 : }
97 :
98 : /// Return the number of frames in each file
99 : /**
100 : * \returns _frames
101 : */
102 0 : int frames()
103 : {
104 0 : return _frames;
105 : }
106 :
107 0 : realT wfPS()
108 : {
109 0 : return _wfPS;
110 : }
111 :
112 0 : void wfPS( realT ps )
113 : {
114 0 : _wfPS = ps;
115 0 : calcPixVal();
116 0 : }
117 :
118 0 : realT F0Photons()
119 : {
120 0 : return _F0Photons;
121 : }
122 :
123 0 : void F0Photons( realT f0 )
124 : {
125 0 : _F0Photons = f0;
126 0 : calcPixVal();
127 0 : }
128 :
129 0 : realT starMag()
130 : {
131 0 : return _starMag;
132 : }
133 :
134 0 : void starMag( realT sm )
135 : {
136 0 : _starMag = sm;
137 0 : calcPixVal();
138 0 : }
139 :
140 0 : realT pixVal()
141 : {
142 0 : return _pixVal;
143 : }
144 :
145 0 : void calcPixVal()
146 : {
147 : // This is the square root of photons/pixel
148 0 : _pixVal = sqrt( _F0Photons ) * pow( 10., -0.2 * _starMag ) * _wfPS; // pow(_wfPS,2);
149 0 : }
150 :
151 : /// Get the file names of the sequence assuming they are stored in a directory.
152 : /**
153 : * \param [in] dir is the directory to search for phase and amplitude files.
154 : * \param [in] max [optional] specifies the maximum number of files to include the file name list.
155 : *
156 : * \retval 0 on success.
157 : * \retval -1 on error.
158 : */
159 : int turbFnames( std::string dir, int max = 0 );
160 :
161 : void openPhaseFrame( int fn );
162 :
163 : void nextPhase( wavefront<realT> &wf );
164 :
165 : void nextWF( wavefront<realT> &wf );
166 : };
167 :
168 : template <typename realT>
169 0 : int turbSequence<realT>::turbFnames( std::string dir, int max )
170 : {
171 0 : ioutils::getFileNames( _phaseFnames, dir, "", ".pha", ".fits" );
172 0 : ioutils::getFileNames( _ampFnames, dir, "", ".amp", ".fits" );
173 :
174 0 : if( _phaseFnames.size() == 0 )
175 : {
176 0 : mxError( "turbSequence", MXE_FILENOTFOUND, "No turbulent phase files found." );
177 0 : return -1;
178 : }
179 :
180 0 : if( _ampFnames.size() == 0 )
181 : {
182 0 : std::cerr << "turbSequence: no turbulent amplitude files found.\n";
183 :
184 0 : _phaseOnly = true;
185 : }
186 :
187 0 : if( max > 0 )
188 : {
189 0 : const size_t maxFiles = static_cast<size_t>( max );
190 :
191 0 : if( _phaseFnames.size() > maxFiles )
192 0 : _phaseFnames.erase( _phaseFnames.begin() + maxFiles, _phaseFnames.end() );
193 :
194 0 : if( _phaseOnly == false && _ampFnames.size() > maxFiles )
195 0 : _ampFnames.erase( _ampFnames.begin() + maxFiles, _ampFnames.end() );
196 : }
197 :
198 0 : _files = _phaseFnames.size();
199 :
200 0 : openPhaseFrame( 0 );
201 :
202 0 : _size = _currPhase.rows();
203 0 : _nPerCube = _currPhase.planes();
204 :
205 0 : _frames = _nPerCube * _phaseFnames.size();
206 :
207 0 : return 0;
208 : }
209 :
210 : template <typename realT>
211 0 : void turbSequence<realT>::openPhaseFrame( int fn )
212 : {
213 0 : fits::fitsFile<float> ff;
214 :
215 0 : std::cout << _phaseFnames[fn] << "\n";
216 :
217 0 : ff.read( _currPhase, _phaseFnames[fn] );
218 0 : ff.close();
219 :
220 0 : if( _phaseOnly == false )
221 : {
222 0 : ff.read( _currAmp, _ampFnames[fn] );
223 0 : ff.close();
224 : }
225 :
226 0 : _currFileNo = fn;
227 0 : _currFrameNo = 0;
228 0 : }
229 :
230 : template <typename realT>
231 0 : void turbSequence<realT>::nextPhase( wavefront<realT> &wf )
232 : {
233 :
234 0 : int Npix = _pupil->sum();
235 :
236 0 : if( _currFrameNo > _nPerCube - 1 )
237 : {
238 0 : openPhaseFrame( _currFileNo + 1 );
239 : }
240 :
241 0 : wf.phase.resize( _pupil->rows(), _pupil->cols() );
242 : // mx::cutImage(wf.phase, _currPhase.image(_currFrameNo), _pupil->rows(), _pupil->cols());
243 :
244 0 : wf.phase = _currPhase.image( _currFrameNo )
245 0 : .block( 0.5 * ( _currPhase.rows() - _pupil->rows() ),
246 0 : 0.5 * ( _currPhase.cols() - _pupil->cols() ),
247 0 : _pupil->rows(),
248 0 : _pupil->cols() );
249 :
250 0 : wf.phase = ( wf.phase - ( wf.phase * ( *_pupil ) ).sum() / Npix ) * ( *_pupil );
251 :
252 : // wf.amplitude.Constant(_pixVal);
253 :
254 0 : wf.amplitude = ( *_pupil ) * _pixVal;
255 :
256 0 : ++_currFrameNo;
257 0 : }
258 :
259 : template <typename realT>
260 0 : void turbSequence<realT>::nextWF( wavefront<realT> &wf )
261 : {
262 :
263 0 : int Npix = _pupil->sum();
264 :
265 0 : if( _currFrameNo > _nPerCube - 1 )
266 : {
267 0 : openPhaseFrame( _currFileNo + 1 );
268 : }
269 0 : wf.phase.resize( _pupil->rows(), _pupil->cols() );
270 : // mx::cutImage(wf.phase, _currPhase.image(_currFrameNo), _pupil->rows(), _pupil->cols());
271 0 : wf.phase = _currPhase.image( _currFrameNo )
272 0 : .block( 0.5 * ( _currPhase.rows() - _pupil->rows() ),
273 0 : 0.5 * ( _currPhase.cols() - _pupil->cols() ),
274 0 : _pupil->rows(),
275 0 : _pupil->cols() );
276 :
277 0 : wf.phase = ( wf.phase - ( wf.phase * ( *_pupil ) ).sum() / Npix ) * ( *_pupil );
278 :
279 0 : if( _phaseOnly )
280 : {
281 0 : wf.amplitude = _pixVal * ( *_pupil );
282 : }
283 : else
284 : {
285 0 : wf.amplitude.resize( _pupil->rows(), _pupil->cols() );
286 : // mx::cutImage(wf.amplitude, _currAmp.image(_currFrameNo), _pupil->rows(), _pupil->cols());
287 0 : wf.amplitude = _currAmp.image( _currFrameNo )
288 0 : .block( 0.5 * ( _currAmp.rows() - _pupil->rows() ),
289 0 : 0.5 * ( _currAmp.cols() - _pupil->cols() ),
290 0 : _pupil->rows(),
291 0 : _pupil->cols() );
292 :
293 0 : wf.amplitude *= _pixVal * ( *_pupil );
294 : }
295 0 : _currFrameNo++;
296 0 : }
297 :
298 : } // namespace sim
299 : } // namespace AO
300 : } // namespace mx
301 :
302 : #endif // turbSequence_hpp
|
