Line data Source code
1 : /** \file idealCoronagraph.hpp
2 : * \brief Declares and defines a class to describe the Ideal Coronagraph.
3 : * \ingroup imaging_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 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 __idealCoronagraph_hpp__
28 : #define __idealCoronagraph_hpp__
29 :
30 : #include <iostream>
31 :
32 : #include "fraunhoferPropagator.hpp"
33 : #include "imagingUtils.hpp"
34 :
35 : #include "../wfp/imagingArray.hpp"
36 :
37 : #include "../ioutils/fits/fitsFile.hpp"
38 : #include "../ioutils/fits/fitsUtils.hpp"
39 : #include "../error/mxErrorOld.hpp"
40 : #include "../improc/imagePads.hpp"
41 :
42 : namespace mx
43 : {
44 :
45 : namespace wfp
46 : {
47 :
48 : /// The Ideal Coronagraph
49 : /** A simple toy coronagraph which operates only the pupil plane, subtracting the energy-minimizing wavefront described
50 : * by Cavarroc et al. \cite cavarroc_2006
51 : * For an on-axis source with a perfectly flat wavefront this results in perfect extinction. This
52 : * coronagraph is not physically realizable, but it is often useful for modeling and analysis, and it has been shown
53 : * that real-world optimized coronagraphs often approach ideal performance \cite sauvage_2010 [and Males in prep].
54 : *
55 : * \ingroup coronagraphs
56 : */
57 : template <typename _realT>
58 : struct idealCoronagraph
59 : {
60 : /// The real floating point type
61 : typedef _realT realT;
62 :
63 : /// The complex floating point type
64 : typedef std::complex<realT> complexT;
65 :
66 : /// The wavefront complex field type
67 : typedef imagingArray<std::complex<realT>, fftwAllocator<std::complex<realT>>, 0> complexFieldT;
68 :
69 : /// The image type
70 : typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
71 :
72 : /// The directory where coronagraph files are stored
73 : std::string _fileDir;
74 :
75 : /// The linear size of the wavefront in pixels
76 : int _wfSz;
77 :
78 : imageT _realPupil;
79 :
80 : complexFieldT m_focalPlane;
81 :
82 : /// Fraunhofer propagator
83 : fraunhoferPropagator<complexFieldT> m_fi;
84 :
85 : /// Determines how the image is centered.
86 : /** If 0 (default) it is at 0.5*(wfSz-1), if true it is shifted by 0.5*m_wholePixel in each axis. This is passed to
87 : * fraunhoferPropagator when it is resized.
88 : */
89 : realT m_wholePixel{ 0 };
90 :
91 : idealCoronagraph();
92 :
93 : /// Get the wavefront size in pixels
94 : /**
95 : * \returns the wavefront size in pixels
96 : */
97 : int wfSz();
98 :
99 : /// Set the wavefront size in pixels.
100 : /**
101 : */
102 : void wfSz( int sz /**< [in] is the new size */ );
103 :
104 : /// Set the real pupil mask.
105 : /** The input mask does not have to be the same size as wfSz, as the stored mask will be padded.
106 : *
107 : * \returns 0 on success.
108 : * \returns -1 on error.
109 : */
110 : int setPupil( imageT &pupil /**< [in] the pupil mask. */ );
111 :
112 : /// Load the real pupil mask from a FITS file.
113 : /** The input mask does not have to be the same size as wfSz, as the stored mask will be padded.
114 : *
115 : * \returns 0 on success.
116 : * \returns -1 on error.
117 : */
118 : int loadPupil( const std::string &pupilFile /**< [in] the FITS file containing the pupil mask. */ );
119 :
120 : /// Load the components of the coronagraph (just a pupil) based in its base name
121 : /** Looks in _filDir for the files.
122 : *
123 : * \returns 0 on success.
124 : * \returns -1 on error.
125 : */
126 : int loadCoronagraph(
127 : const std::string &cName /**< The name of the coronagraph, without directory or file extensions */ );
128 :
129 : /// Propagate the given pupil-plane wavefront through the coronagraph to the exit pupil plane
130 : int propagate(
131 : complexFieldT
132 : &pupilPlane /**< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph. */ );
133 :
134 : /// Propagate the given pupil-plane wavefront through the coronagraph to the exit pupil plane, and then to the final
135 : /// focal plane.
136 : int
137 : propagate( imageT &fpIntensity, ///< [out] The intensity image in the focal plane. This should be pre-allocated.
138 : complexFieldT
139 : &pupilPlane ///< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph.
140 : );
141 :
142 : /// Propagate the given pupil-plane wavefront without the coronagraph.
143 : /** For the ideal coronagraph nothing is done. This method is included for compliance with
144 : * with the coronagraph interface.
145 : */
146 : int propagateNC(
147 : complexFieldT &pupilPlane /**< [in.out] The wavefront at the input pupil plane. It is un-modified. */ );
148 :
149 : /// Propagate the given pupil-plane wavefront without the coronagraph to the exit pupil plane, and then to the final
150 : /// focal plane.
151 : /** For the ideal coronagraph nothing is done to the input wavefront.
152 : */
153 : int
154 : propagateNC( imageT &fpIntensity, ///< [out] The intensity image in the focal plane. This should be pre-allocated.
155 : complexFieldT &pupilPlane ///< [in.out] The wavefront at the input pupil plane. It is un-modified.
156 : );
157 :
158 : bool apodize{ false };
159 : imageT apodizer;
160 : };
161 :
162 : template <typename realT>
163 0 : idealCoronagraph<realT>::idealCoronagraph()
164 : {
165 0 : _wfSz = 0;
166 0 : }
167 :
168 : template <typename realT>
169 0 : int idealCoronagraph<realT>::wfSz()
170 : {
171 0 : return _wfSz;
172 : }
173 :
174 : template <typename realT>
175 0 : void idealCoronagraph<realT>::wfSz( int sz )
176 : {
177 0 : _wfSz = sz;
178 :
179 0 : m_fi.setWavefrontSizePixels( sz );
180 0 : m_fi.wholePixel( m_wholePixel );
181 0 : m_focalPlane.resize( sz, sz );
182 0 : }
183 :
184 : template <typename realT>
185 0 : int idealCoronagraph<realT>::setPupil( imageT &pupil )
186 : {
187 0 : if( _wfSz <= 0 )
188 : {
189 0 : mxError( "idealCoronagraph", MXE_PARAMNOTSET, "Must set wavefront size (wfSz) before setting up coronagraph." );
190 0 : return -1;
191 : }
192 :
193 : // Create Coronagraph pupil.
194 0 : improc::padImage( _realPupil, pupil, 0.5 * ( _wfSz - pupil.rows() ), 0 );
195 :
196 0 : return 0;
197 : }
198 :
199 : template <typename realT>
200 0 : int idealCoronagraph<realT>::loadPupil( const std::string &pupilFile )
201 : {
202 :
203 0 : imageT pupil;
204 :
205 0 : fits::fitsFile<realT> ff;
206 :
207 0 : ff.read( pupil, pupilFile );
208 :
209 0 : return setPupil( pupil );
210 :
211 : return 0;
212 0 : }
213 :
214 : template <typename realT>
215 0 : int idealCoronagraph<realT>::loadCoronagraph( const std::string &cName )
216 : {
217 :
218 0 : if( _fileDir == "" )
219 : {
220 0 : mxError( "idealCoronagraph", MXE_PARAMNOTSET, "file directory (fileDir) not set." );
221 0 : return -1;
222 : }
223 :
224 0 : std::string pupilFile = _fileDir + "/" + cName + ".fits";
225 :
226 0 : return loadPupil( pupilFile );
227 0 : }
228 :
229 : template <typename realT>
230 0 : int idealCoronagraph<realT>::propagate( complexFieldT &pupilPlane )
231 : {
232 0 : if( pupilPlane.rows() != _realPupil.rows() || pupilPlane.cols() != _realPupil.cols() )
233 : {
234 0 : mxError( "idealCoronagraph", MXE_SIZEERR, "pupilPlane wavefront size does not match realPupil" );
235 0 : return -1;
236 : }
237 :
238 0 : Eigen::Map<Eigen::Array<complexT, -1, -1>> eigWf( pupilPlane.data(), pupilPlane.cols(), pupilPlane.rows() );
239 0 : Eigen::Map<Eigen::Array<realT, -1, -1>> eigPup( _realPupil.data(), _realPupil.cols(), _realPupil.rows() );
240 :
241 0 : eigWf = ( eigWf - ( ( eigWf * eigPup ).sum() / ( eigPup * eigPup ).sum() ) ) * eigPup;
242 :
243 0 : if( apodize )
244 : {
245 0 : eigWf *= apodizer;
246 : }
247 :
248 0 : return 0;
249 : }
250 :
251 : template <typename realT>
252 0 : int idealCoronagraph<realT>::propagate( imageT &fpIntensity, complexFieldT &pupilPlane )
253 : {
254 0 : propagate( pupilPlane );
255 :
256 0 : m_fi.propagatePupilToFocal( m_focalPlane, pupilPlane );
257 :
258 0 : int x0 = 0.5 * ( _wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
259 0 : int y0 = 0.5 * ( _wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
260 :
261 0 : extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
262 :
263 0 : return 0;
264 : }
265 :
266 : template <typename realT>
267 0 : int idealCoronagraph<realT>::propagateNC( complexFieldT &pupilPlane )
268 : {
269 : static_cast<void>( pupilPlane );
270 0 : return 0;
271 : }
272 :
273 : template <typename realT>
274 0 : int idealCoronagraph<realT>::propagateNC( imageT &fpIntensity, complexFieldT &pupilPlane )
275 : {
276 0 : propagateNC( pupilPlane );
277 :
278 0 : m_fi.propagatePupilToFocal( m_focalPlane, pupilPlane );
279 :
280 0 : int x0 = 0.5 * ( _wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
281 0 : int y0 = 0.5 * ( _wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
282 :
283 0 : extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
284 :
285 0 : return 0;
286 : }
287 :
288 : } // namespace wfp
289 : } // namespace mx
290 :
291 : #endif //__idealCoronagraph_hpp__
|
