Line data Source code
1 : /** \file lyotCoronagraph.hpp
2 : * \brief Declares and defines a class to describe and optimize a Lyot Coronagraph.
3 : * \ingroup imaging_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017, 2018 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 __lyotCoronagraph_hpp__
28 : #define __lyotCoronagraph_hpp__
29 :
30 : #include <iostream>
31 :
32 : #include "imagingArray.hpp"
33 : #include "imagingUtils.hpp"
34 : #include "fraunhoferPropagator.hpp"
35 :
36 : #include "../ioutils/fits/fitsFile.hpp"
37 : #include "../ioutils/fits/fitsUtils.hpp"
38 :
39 : #include "../error/mxErrorOld.hpp"
40 :
41 : #include "../improc/eigenImage.hpp"
42 : #include "../improc/eigenCube.hpp"
43 :
44 : namespace mx
45 : {
46 :
47 : namespace wfp
48 : {
49 :
50 : /// The Lyot Coronagraph
51 : /** A generalized Lyot coronagraph, which can include a pupil apodization, a focal plane mask with complex transmission,
52 : * and a Lyot stop. Light is propagated through the coronagraph plane-by-plane, and the complex wavefront can be
53 : * accessed at any plane. Also provides functions for optimizing the pupil apodizer and focal plane mask.
54 : *
55 : * \ingroup coronagraphs
56 : */
57 : template <typename _realT, typename _fpmaskFloatT>
58 : struct lyotCoronagraph
59 : {
60 : public:
61 : /// The real floating point type
62 : typedef _realT realT;
63 :
64 : /// The real floating point type for mask calculations.
65 : typedef _fpmaskFloatT fpmaskFloatT;
66 :
67 : /// The complex floating point type
68 : typedef std::complex<realT> complexT;
69 :
70 : /// The wavefront complex field type
71 : typedef mx::wfp::imagingArray<std::complex<realT>, fftwAllocator<std::complex<realT>>, 0> complexFieldT;
72 :
73 : /// The focal plane mask type
74 : typedef Eigen::Array<std::complex<fpmaskFloatT>, Eigen::Dynamic, Eigen::Dynamic> fpMaskT;
75 :
76 : /// The image type
77 : typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
78 :
79 : public:
80 : /// The directory where coronagraph files are stored
81 : std::string m_fileDir{ "coron" };
82 :
83 : /// The linear size of the wavefront in pixels
84 : int m_wfSz{ 0 };
85 :
86 : /// Image containing the pupil apodization.
87 : imageT m_pupilApodizer;
88 :
89 : int m_maskSource; ///< 0= read from file, 1 = constructed by makeFocalMask, 2 = trans optimized.
90 : std::string m_maskFile; ///< Name of file from which mask was loaded.
91 :
92 : realT m_maskRad; ///< Radius of mask if it was constructed.
93 : realT m_maskTrans; ///< Transmission of mask if it was constructed.
94 :
95 : /// The focal plane mask.
96 : fpMaskT m_focalMask;
97 :
98 : /// Image containing the lyot stop.
99 : imageT m_lyotStop;
100 :
101 : complexFieldT m_focalPlane;
102 :
103 : bool m_savePreMaskFocalPlane{ false };
104 : complexFieldT m_preMaskFocalPlane;
105 :
106 : bool m_savePreLyotPupilPlane{ false };
107 : complexFieldT m_preLyotPupilPlane;
108 :
109 : /// Fraunhofer propagator
110 : fraunhoferPropagator<complexFieldT> m_fp;
111 :
112 : public:
113 : /// Default c'tor
114 : lyotCoronagraph();
115 :
116 : /// Get the wavefront size in pixels
117 : /**
118 : * \returns the wavefront size in pixels
119 : */
120 : int wfSz();
121 :
122 : /// Set the wavefront size in pixels.
123 : /**
124 : */
125 : void wfSz( int sz /**< [in] is the new size */ );
126 :
127 : /// Make the focal plane mask
128 : void makeFocalMask( realT rad, fpmaskFloatT trans = 0.0, int sz = 0.0 );
129 :
130 : /// Load the apodizer from a FITS file
131 : /**
132 : * \returns 0 on success.
133 : * \returns -1 on error.
134 : */
135 : int loadApodizer( const std::string &apodName /**< [in] is the name of the FITS file containing the apodizer. */ );
136 :
137 : /// Load the focal plane mask from a FITS file
138 : /**
139 : * \returns 0 on success.
140 : * \returns -1 on error.
141 : */
142 : int loadFocalMask(
143 : const std::string &fpmName /**< [in] is the name of the FITS file containing the focal plane mask. */ );
144 :
145 : /// Load the Lyot stop from a FITS file
146 : /**
147 : * \returns 0 on success.
148 : * \returns -1 on error.
149 : */
150 : int loadLyotStop( const std::string &lyotName /**< [in] is the name of the FITS file containing the lyot stop. */ );
151 :
152 : /// Load the components of the coronagraph from FITS files
153 : /**
154 : * \returns 0 on success.
155 : * \returns -1 on error.
156 : */
157 : int
158 : loadCoronagraph( const std::string &apodName, ///< [in] is the name of the FITS file containing the apodizer.
159 : const std::string &fpmName, ///< [in] is the name of the FITS file containing the focal plane mask.
160 : const std::string &lyotName ///< [in] is the name of the FITS file containing the lyot stop.
161 : );
162 :
163 : /// Load the components of the coronagraph based in its base name
164 : /** Looks in _filDir for the files.
165 : *
166 : * \returns 0 on success.
167 : * \returns -1 on error.
168 : */
169 : int
170 : loadCoronagraph( const std::string
171 : &cName /**< [in] is the base name of the coronagraph without directory or file extension. */ );
172 :
173 : void applyApodizer( complexFieldT &pupilPlane );
174 : void applyFocalMask( complexFieldT &focalPlane );
175 : void applyLyotStop( complexFieldT &lyotPlane );
176 :
177 : /// Propagate the given pupil-plane wavefront through the coronagraph.
178 : int propagate(
179 : complexFieldT
180 : &pupilPlane /**< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph. */ );
181 :
182 : int
183 : propagate( imageT &fpIntensity, ///< [out] The intensity image in the focal plane. This should be pre-allocated.
184 : complexFieldT
185 : &pupilPlane ///< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph.
186 : );
187 :
188 : /// Propagate the given pupil-plane wavefront without the coronagraph.
189 : /** For a Lyot coronagraph, this applies the pupil apodization and Lyot stop, but not the FPM, to the given
190 : * pupil-plane wavefront such that the result will produce the non-coronagraphic (off-axis) PSF.
191 : */
192 : int propagateNC(
193 : complexFieldT
194 : &pupilPlane /**< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph. */ );
195 :
196 : int propagateNC(
197 : imageT &fpIntensity, ///< [out] The intensity image in the focal plane. This should be pre-allocated.
198 : complexFieldT
199 : &pupilPlane ///< [in.out] The wavefront at the input pupil plane. It is modified by the coronagraph.
200 : );
201 :
202 : /// Optimize the pupil amplitude apodization and focal-plane mask complex transmission.
203 : /** Uses the algorithm in Guyon (2014) \cite guyon_2014 for optimizing an Apodized-Pupil Lyot Complex Mask
204 : * Coronagraph (APLCMC).
205 : *
206 : */
207 : void
208 : optimizeAPLCMC( imageT &geomPupil, ///< The geometric pupil mask, binary 1/0 transmission.
209 : realT fpmRadPix, ///< The radius in pixels of the FPM.
210 : realT relTol, ///< Relative tolerance for convergence.
211 : realT absTol, ///< Absolute tolerance for convergence.
212 : int maxIter, ///< Maximum number of iterations to allow.
213 : const std::string &cname ///< Name of coronagraph, used as base-name for output files in m_fileDir.
214 : );
215 :
216 : void
217 : optimizeAPLCMC( imageT &geomPupil, ///< The geometric pupil mask, binary 1/0 transmission.
218 : realT fpmRadPix, ///< The radius in pixels of the FPM.
219 : imageT &fpmPhase,
220 : realT relTol, ///< Relative tolerance for convergence.
221 : realT absTol, ///< Absolute tolerance for convergence.
222 : int maxIter, ///< Maximum number of iterations to allow.
223 : const std::string &cname ///< Name of coronagraph, used as base-name for output files in m_fileDir.
224 : );
225 : };
226 :
227 : template <typename _realT, typename _fpmaskFloatT>
228 0 : lyotCoronagraph<_realT, _fpmaskFloatT>::lyotCoronagraph()
229 : {
230 0 : }
231 :
232 : template <typename _realT, typename _fpmaskFloatT>
233 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::wfSz()
234 : {
235 0 : return m_wfSz;
236 : }
237 :
238 : template <typename _realT, typename _fpmaskFloatT>
239 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::wfSz( int sz )
240 : {
241 0 : m_wfSz = sz;
242 :
243 0 : m_fp.setWavefrontSizePixels( sz );
244 0 : m_focalPlane.resize( sz, sz );
245 0 : }
246 :
247 : template <typename _realT, typename _fpmaskFloatT>
248 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::makeFocalMask( _realT rad, _fpmaskFloatT trans, int sz )
249 : {
250 0 : if( sz == 0 )
251 : {
252 0 : sz = 2 * rad + 1;
253 0 : if( sz % 2 == 1 )
254 0 : ++sz;
255 : }
256 :
257 : // Make a real circular mask
258 0 : improc::eigenImage<fpmaskFloatT> fpm( sz, sz );
259 0 : mx::wfp::circularPupil( fpm, 0., rad );
260 :
261 : // Convert to complex amplitude
262 0 : m_focalMask.resize( sz, sz );
263 0 : for( int i = 0; i < sz; ++i )
264 : {
265 0 : for( int j = 0; j < sz; ++j )
266 : {
267 0 : if( fpm( i, j ) == 1 )
268 0 : m_focalMask( i, j ) = std::complex<fpmaskFloatT>( trans, 0 );
269 : else
270 0 : m_focalMask( i, j ) = std::complex<fpmaskFloatT>( 1, 0 );
271 : }
272 : }
273 :
274 0 : m_maskSource = 1;
275 0 : m_maskFile = "";
276 0 : m_maskRad = rad;
277 0 : m_maskTrans = 0.0;
278 0 : }
279 :
280 : template <typename _realT, typename _fpmaskFloatT>
281 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::loadApodizer( const std::string &apodName )
282 : {
283 0 : fits::fitsFile<_realT> ff;
284 :
285 0 : if( !!ff.read( m_pupilApodizer, apodName ) )
286 0 : return -1;
287 :
288 0 : return 0;
289 0 : }
290 :
291 : template <typename _realT, typename _fpmaskFloatT>
292 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::loadFocalMask( const std::string &fpmName )
293 : {
294 0 : fits::fitsFile<_realT> ff;
295 0 : improc::eigenCube<fpmaskFloatT> fpm;
296 :
297 0 : if( !!ff.read( fpm, fpmName ) )
298 0 : return -1;
299 :
300 0 : if( fpm.planes() == 1 )
301 : {
302 0 : m_focalMask = fpm.image( 0 );
303 : }
304 0 : else if( fpm.planes() == 2 )
305 : {
306 0 : m_focalMask.resize( fpm.rows(), fpm.cols() );
307 :
308 0 : for( int r = 0; r < fpm.rows(); ++r )
309 : {
310 0 : for( int c = 0; c < fpm.cols(); ++c )
311 : {
312 0 : m_focalMask( r, c ) =
313 0 : fpm.image( 0 )( r, c ) * exp( std::complex<fpmaskFloatT>( 0, fpm.image( 1 )( r, c ) ) );
314 : }
315 : }
316 : }
317 : else
318 : {
319 0 : std::cerr << "too many planes in focal mask file\n";
320 0 : return -1;
321 : }
322 :
323 0 : m_maskSource = 0;
324 0 : m_maskFile = fpmName;
325 0 : m_maskRad = 0.0;
326 0 : m_maskTrans = 0.0;
327 :
328 0 : return 0;
329 0 : }
330 :
331 : template <typename _realT, typename _fpmaskFloatT>
332 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::loadLyotStop( const std::string &lyotName )
333 : {
334 0 : fits::fitsFile<_realT> ff;
335 :
336 0 : if( !!ff.read( m_lyotStop, lyotName ) )
337 0 : return -1;
338 :
339 0 : return 0;
340 0 : }
341 :
342 : template <typename _realT, typename _fpmaskFloatT>
343 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::loadCoronagraph( const std::string &apodName,
344 : const std::string &fpmName,
345 : const std::string &lyotName )
346 : {
347 0 : if( loadApodizer( apodName ) < 0 )
348 0 : return -1;
349 0 : if( loadFocalMask( fpmName ) < 0 )
350 0 : return -1;
351 0 : if( loadLyotStop( lyotName ) < 0 )
352 0 : return -1;
353 :
354 0 : return 0;
355 : }
356 :
357 : template <typename _realT, typename _fpmaskFloatT>
358 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::loadCoronagraph( const std::string &cName )
359 : {
360 0 : if( m_fileDir == "" )
361 : {
362 0 : mxError( "lyotCoronagraph", MXE_PARAMNOTSET, "file directory (fileDir) not set." );
363 0 : return -1;
364 : }
365 :
366 0 : std::string apodName = m_fileDir + "/" + cName + "_apod.fits";
367 0 : std::string fpmName = m_fileDir + "/" + cName + "_fpm.fits";
368 0 : std::string lyotName = m_fileDir + "/" + cName + "_lyot.fits";
369 :
370 0 : return loadCoronagraph( apodName, fpmName, lyotName );
371 0 : }
372 :
373 : template <typename _realT, typename _fpmaskFloatT>
374 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::applyApodizer( complexFieldT &pupilPlane )
375 : {
376 0 : int sz = m_pupilApodizer.rows();
377 :
378 0 : realT w = 0.5 * ( sz - 1.0 );
379 :
380 0 : realT xc = 0.5 * ( pupilPlane.rows() - 1 );
381 0 : realT yc = 0.5 * ( pupilPlane.cols() - 1 );
382 :
383 0 : for( int i = 0; i < pupilPlane.rows(); ++i )
384 : {
385 0 : for( int j = 0; j < pupilPlane.cols(); ++j )
386 : {
387 0 : if( i >= xc - w && i <= xc + w && j >= yc - w && j <= yc + w )
388 0 : pupilPlane( i, j ) *= m_pupilApodizer( (int)( i - ( xc - w ) ), (int)( j - ( yc - w ) ) );
389 : else
390 0 : pupilPlane( i, j ) *= 0;
391 : }
392 : }
393 0 : }
394 :
395 : template <typename _realT, typename _fpmaskFloatT>
396 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::applyFocalMask( complexFieldT &focalPlane )
397 : {
398 0 : int sz = m_focalMask.rows();
399 :
400 0 : realT w = 0.5 * ( sz - 1.0 );
401 :
402 0 : realT xc = 0.5 * ( focalPlane.rows() - 1 );
403 0 : realT yc = 0.5 * ( focalPlane.cols() - 1 );
404 :
405 0 : for( int i = 0; i < sz; ++i )
406 : {
407 0 : for( int j = 0; j < sz; ++j )
408 : {
409 0 : focalPlane( xc - w + i, yc - w + j ) *= m_focalMask( i, j );
410 : }
411 : }
412 0 : }
413 :
414 : template <typename _realT, typename _fpmaskFloatT>
415 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::applyLyotStop( complexFieldT &lyotPlane )
416 : {
417 0 : int sz = m_lyotStop.rows();
418 :
419 0 : realT w = 0.5 * ( sz - 1.0 );
420 :
421 0 : realT xc = 0.5 * ( lyotPlane.rows() - 1 );
422 0 : realT yc = 0.5 * ( lyotPlane.cols() - 1 );
423 :
424 0 : for( int i = 0; i < lyotPlane.rows(); ++i )
425 : {
426 0 : for( int j = 0; j < lyotPlane.cols(); ++j )
427 : {
428 0 : if( i >= xc - w && i <= xc + w && j >= yc - w && j <= yc + w )
429 0 : lyotPlane( i, j ) *= m_lyotStop( (int)( i - ( xc - w ) ), (int)( j - ( yc - w ) ) );
430 : else
431 0 : lyotPlane( i, j ) = 0;
432 : }
433 : }
434 0 : }
435 :
436 : template <typename _realT, typename _fpmaskFloatT>
437 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::propagate( complexFieldT &pupilPlane )
438 : {
439 0 : applyApodizer( pupilPlane );
440 :
441 0 : m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
442 :
443 0 : if( m_savePreMaskFocalPlane )
444 : {
445 0 : m_preMaskFocalPlane = m_focalPlane;
446 : }
447 :
448 0 : applyFocalMask( m_focalPlane );
449 :
450 0 : m_fp.propagateFocalToPupil( pupilPlane, m_focalPlane );
451 :
452 0 : if( m_savePreLyotPupilPlane )
453 : {
454 0 : m_preLyotPupilPlane = pupilPlane;
455 : }
456 :
457 0 : applyLyotStop( pupilPlane );
458 :
459 0 : return 0;
460 : }
461 :
462 : template <typename _realT, typename _fpmaskFloatT>
463 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::propagate( imageT &fpIntensity, complexFieldT &pupilPlane )
464 : {
465 0 : propagate( pupilPlane );
466 :
467 0 : m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
468 :
469 0 : int x0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
470 0 : int y0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
471 :
472 0 : extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
473 :
474 0 : return 0;
475 : }
476 : template <typename _realT, typename _fpmaskFloatT>
477 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::propagateNC( complexFieldT &pupilPlane )
478 : {
479 0 : applyApodizer( pupilPlane );
480 :
481 0 : m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
482 :
483 0 : m_fp.propagateFocalToPupil( pupilPlane, m_focalPlane );
484 :
485 0 : applyLyotStop( pupilPlane );
486 :
487 0 : return 0;
488 : }
489 :
490 : template <typename _realT, typename _fpmaskFloatT>
491 0 : int lyotCoronagraph<_realT, _fpmaskFloatT>::propagateNC( imageT &fpIntensity, complexFieldT &pupilPlane )
492 : {
493 0 : propagateNC( pupilPlane );
494 :
495 0 : m_fp.propagatePupilToFocal( m_focalPlane, pupilPlane );
496 :
497 0 : int x0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.rows() - 1 );
498 0 : int y0 = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( fpIntensity.cols() - 1 );
499 :
500 0 : extractIntensityImage( fpIntensity, 0, fpIntensity.rows(), 0, fpIntensity.cols(), m_focalPlane, x0, y0 );
501 :
502 0 : return 0;
503 : }
504 :
505 : template <typename _realT, typename _fpmaskFloatT>
506 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::optimizeAPLCMC(
507 : imageT &geomPupil, realT fpmRadPix, realT relTol, realT absTol, int maxIter, const std::string &cname )
508 : {
509 0 : complexFieldT focalPlane, pupilPlane;
510 :
511 0 : pupilPlane.resize( m_wfSz, m_wfSz );
512 0 : focalPlane.resize( m_wfSz, m_wfSz );
513 :
514 0 : mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> mask( m_wfSz, m_wfSz );
515 0 : mx::wfp::circularPupil( mask, 0., fpmRadPix );
516 : /*for(int c=0; c< mask.cols()*0.5; ++c)
517 : {
518 : for(int r=0; r < mask.rows(); ++r) mask(r,c) = 0;
519 : }*/
520 :
521 : // Initialize pupilImage
522 0 : mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> pupilImage( m_wfSz, m_wfSz );
523 0 : pupilImage.setZero();
524 :
525 0 : int gpLLi = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.rows() - 1 );
526 0 : int gpLLj = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.cols() - 1 );
527 :
528 0 : int gpURi = gpLLi + geomPupil.rows();
529 0 : int gpURj = gpLLj + geomPupil.cols();
530 :
531 0 : for( int i = 0; i < geomPupil.rows(); ++i )
532 : {
533 0 : for( int j = 0; j < geomPupil.cols(); ++j )
534 : {
535 0 : pupilImage( gpLLi + i, gpLLj + j ) = geomPupil( i, j );
536 : }
537 : }
538 :
539 : realT lastLambdaA, LambdaA;
540 :
541 0 : lastLambdaA = 1;
542 : int n;
543 0 : std::string reason;
544 0 : for( n = 0; n < maxIter; ++n )
545 : {
546 0 : mx::wfp::makeComplexPupil( pupilPlane, pupilImage, m_wfSz );
547 :
548 0 : m_fp.propagatePupilToFocal( focalPlane, pupilPlane );
549 :
550 0 : for( int i = 0; i < m_wfSz; ++i )
551 : {
552 0 : for( int j = 0; j < m_wfSz; ++j )
553 : {
554 0 : focalPlane( i, j ) *= mask( i, j );
555 : }
556 : }
557 :
558 0 : m_fp.propagateFocalToPupil( pupilPlane, focalPlane );
559 :
560 0 : for( int i = 0; i < m_wfSz; ++i )
561 0 : for( int j = 0; j < m_wfSz; ++j )
562 0 : pupilImage( i, j ) = abs( pupilPlane( i, j ) );
563 :
564 0 : LambdaA = 0;
565 0 : for( int i = 0; i < m_wfSz; ++i )
566 : {
567 0 : for( int j = 0; j < m_wfSz; ++j )
568 : {
569 0 : if( i >= gpLLi && i < gpURi && j >= gpLLj && j < gpURj )
570 : {
571 0 : pupilImage( i, j ) *= geomPupil( i - gpLLi, j - gpLLj );
572 0 : if( pupilImage( i, j ) > LambdaA )
573 0 : LambdaA = pupilImage( i, j );
574 : }
575 : else
576 : {
577 0 : pupilImage( i, j ) = 0;
578 : }
579 : }
580 : }
581 : // LambdaA = 1.0/LambdaA;
582 0 : for( int i = 0; i < m_wfSz; ++i )
583 0 : for( int j = 0; j < m_wfSz; ++j )
584 0 : pupilImage( i, j ) /= LambdaA;
585 :
586 0 : std::cout << n << " " << LambdaA << "\n";
587 0 : if( fabs( LambdaA - lastLambdaA ) < absTol )
588 : {
589 0 : std::cout << "Converged on absTol.\n";
590 0 : reason = "absTol";
591 0 : break;
592 : }
593 0 : if( fabs( ( LambdaA - lastLambdaA ) / lastLambdaA ) < relTol )
594 : {
595 0 : std::cout << "Converged on relTol.\n";
596 0 : reason = "relTol";
597 0 : break;
598 : }
599 :
600 0 : if( n == maxIter - 1 )
601 : {
602 0 : std::cout << "maxIter reached.\n";
603 0 : reason = "maxIter";
604 : }
605 :
606 0 : lastLambdaA = LambdaA;
607 : }
608 0 : if( LambdaA > 1 )
609 0 : LambdaA = 1;
610 :
611 0 : std::cout << "LambdaA: = " << LambdaA << "\n";
612 :
613 0 : realT trans = 1.0 - 1.0 / LambdaA;
614 :
615 0 : int pupSize = geomPupil.rows();
616 :
617 0 : m_pupilApodizer.resize( pupSize, pupSize );
618 :
619 0 : extractBlock( m_pupilApodizer,
620 : 0,
621 : pupSize,
622 : 0,
623 : pupSize,
624 : pupilImage,
625 0 : 0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ),
626 0 : 0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ) );
627 :
628 0 : makeFocalMask( fpmRadPix, trans, pupSize );
629 :
630 : /*for(int c=0; c< m_focalMask.cols()*0.5; ++c)
631 : {
632 : for(int r=0; r < m_focalMask.rows(); ++r)
633 : {
634 : if(m_focalMask(r,c).real()==1) m_focalMask(r,c) = 0;
635 : }
636 : }*/
637 :
638 0 : m_maskSource = 2;
639 :
640 0 : m_lyotStop = geomPupil;
641 :
642 0 : fits::fitsHeader head;
643 :
644 0 : head.append( "", fits::fitsCommentType(), "----------------------------------------" );
645 0 : head.append( "", fits::fitsCommentType(), "lyotCoronagraph optimization Parameters:" );
646 0 : head.append( "", fits::fitsCommentType(), "----------------------------------------" );
647 0 : head.append<int>( "WFSZ", m_wfSz, "Size of wavefront used for FFTs (pixels)" );
648 0 : head.append<realT>( "FPMRADPX", fpmRadPix, "input radius of focal plane mask (pixels)" );
649 0 : head.append<realT>( "ABSTOL", absTol, "input absolute tolerance" );
650 0 : head.append<realT>( "RELTOL", relTol, "input relative tolerance" );
651 0 : head.append<int>( "MAXITER", maxIter, "input maximum iterations" );
652 0 : head.append<int>( "NITER", n, "actual number of iterations" );
653 0 : head.append<std::string>( "XREASON", reason, "reason for convergence" );
654 0 : head.append<realT>( "FPMTRANS", trans, "transmission of FPM" );
655 :
656 0 : fits::fitsFile<double> ff;
657 :
658 0 : std::string fname = "coron/" + cname + "_apod.fits";
659 :
660 0 : ff.write( fname, m_pupilApodizer, head );
661 :
662 0 : fname = "coron/" + cname + "_fpm.fits";
663 0 : improc::eigenImage<fpmaskFloatT> fpm( m_focalMask.rows(), m_focalMask.cols() );
664 0 : for( int r = 0; r < m_focalMask.rows(); ++r )
665 : {
666 0 : for( int c = 0; c < m_focalMask.cols(); ++c )
667 : {
668 0 : fpm( r, c ) = m_focalMask( r, c ).real();
669 : }
670 : }
671 :
672 0 : ff.write( fname, fpm, head );
673 :
674 0 : fname = "coron/" + cname + "_lyot.fits";
675 0 : ff.write( fname, m_lyotStop, head );
676 0 : }
677 :
678 : template <typename _realT, typename _fpmaskFloatT>
679 0 : void lyotCoronagraph<_realT, _fpmaskFloatT>::optimizeAPLCMC( imageT &geomPupil,
680 : realT fpmRadPix,
681 : imageT &fpmPhase,
682 : realT relTol,
683 : realT absTol,
684 : int maxIter,
685 : const std::string &cname )
686 : {
687 0 : complexFieldT focalPlane, pupilPlane;
688 :
689 0 : pupilPlane.resize( m_wfSz, m_wfSz );
690 0 : focalPlane.resize( m_wfSz, m_wfSz );
691 :
692 0 : mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> mask( m_wfSz, m_wfSz );
693 0 : mx::wfp::circularPupil( mask, 0., fpmRadPix );
694 :
695 : // Initialize pupilImage
696 0 : mx::wfp::imagingArray<realT, fftwAllocator<realT>, 0> pupilImage( m_wfSz, m_wfSz );
697 0 : pupilImage.setZero();
698 :
699 0 : int gpLLi = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.rows() - 1 );
700 0 : int gpLLj = 0.5 * ( m_wfSz - 1 ) - 0.5 * ( geomPupil.cols() - 1 );
701 :
702 0 : int gpURi = gpLLi + geomPupil.rows();
703 0 : int gpURj = gpLLj + geomPupil.cols();
704 :
705 0 : for( int i = 0; i < geomPupil.rows(); ++i )
706 : {
707 0 : for( int j = 0; j < geomPupil.cols(); ++j )
708 : {
709 0 : pupilImage( gpLLi + i, gpLLj + j ) = geomPupil( i, j );
710 : }
711 : }
712 :
713 : realT lastLambdaA, LambdaA;
714 :
715 0 : lastLambdaA = 1;
716 : int n;
717 0 : std::string reason;
718 0 : for( n = 0; n < maxIter; ++n )
719 : {
720 0 : mx::wfp::makeComplexPupil( pupilPlane, pupilImage, m_wfSz );
721 :
722 0 : m_fp.propagatePupilToFocal( focalPlane, pupilPlane );
723 :
724 0 : for( int i = 0; i < m_wfSz; ++i )
725 : {
726 0 : for( int j = 0; j < m_wfSz; ++j )
727 : {
728 0 : focalPlane( i, j ) *= mask( i, j ) * exp( std::complex<realT>( 0, fpmPhase( i, j ) ) );
729 : ;
730 : }
731 : }
732 :
733 0 : m_fp.propagateFocalToPupil( pupilPlane, focalPlane );
734 :
735 0 : for( int i = 0; i < m_wfSz; ++i )
736 0 : for( int j = 0; j < m_wfSz; ++j )
737 0 : pupilImage( i, j ) = abs( pupilPlane( i, j ) );
738 :
739 0 : LambdaA = 0;
740 0 : for( int i = 0; i < m_wfSz; ++i )
741 : {
742 0 : for( int j = 0; j < m_wfSz; ++j )
743 : {
744 0 : if( i >= gpLLi && i < gpURi && j >= gpLLj && j < gpURj )
745 : {
746 0 : pupilImage( i, j ) *= geomPupil( i - gpLLi, j - gpLLj );
747 0 : if( pupilImage( i, j ) > LambdaA )
748 0 : LambdaA = pupilImage( i, j );
749 : }
750 : else
751 : {
752 0 : pupilImage( i, j ) = 0;
753 : }
754 : }
755 : }
756 : // LambdaA = 1.0/LambdaA;
757 0 : for( int i = 0; i < m_wfSz; ++i )
758 0 : for( int j = 0; j < m_wfSz; ++j )
759 0 : pupilImage( i, j ) /= LambdaA;
760 :
761 0 : std::cout << n << " " << LambdaA << "\n";
762 0 : if( fabs( LambdaA - lastLambdaA ) < absTol )
763 : {
764 0 : std::cout << "Converged on absTol.\n";
765 0 : reason = "absTol";
766 0 : break;
767 : }
768 0 : if( fabs( ( LambdaA - lastLambdaA ) / lastLambdaA ) < relTol )
769 : {
770 0 : std::cout << "Converged on relTol.\n";
771 0 : reason = "relTol";
772 0 : break;
773 : }
774 :
775 0 : if( n == maxIter - 1 )
776 : {
777 0 : std::cout << "maxIter reached.\n";
778 0 : reason = "maxIter";
779 : }
780 :
781 0 : lastLambdaA = LambdaA;
782 : }
783 0 : if( LambdaA > 1 )
784 0 : LambdaA = 1;
785 :
786 0 : std::cout << "LambdaA: = " << LambdaA << "\n";
787 :
788 0 : realT trans = 1.0 - 1.0 / LambdaA;
789 :
790 0 : int pupSize = geomPupil.rows();
791 :
792 0 : m_pupilApodizer.resize( pupSize, pupSize );
793 :
794 0 : extractBlock( m_pupilApodizer,
795 : 0,
796 : pupSize,
797 : 0,
798 : pupSize,
799 : pupilImage,
800 0 : 0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ),
801 0 : 0.5 * ( ( pupilImage.rows() - 1 ) - ( pupSize - 1 ) ) );
802 :
803 0 : makeFocalMask( fpmRadPix, trans, pupSize );
804 0 : m_maskSource = 2;
805 :
806 0 : m_lyotStop = geomPupil;
807 :
808 0 : fits::fitsHeader head;
809 :
810 0 : head.append( "", fits::fitsCommentType(), "----------------------------------------" );
811 0 : head.append( "", fits::fitsCommentType(), "lyotCoronagraph optimization Parameters:" );
812 0 : head.append( "", fits::fitsCommentType(), "----------------------------------------" );
813 0 : head.append<int>( "WFSZ", m_wfSz, "Size of wavefront used for FFTs (pixels)" );
814 0 : head.append<realT>( "FPMRADPX", fpmRadPix, "input radius of focal plane mask (pixels)" );
815 0 : head.append<realT>( "ABSTOL", absTol, "input absolute tolerance" );
816 0 : head.append<realT>( "RELTOL", relTol, "input relative tolerance" );
817 0 : head.append<int>( "MAXITER", maxIter, "input maximum iterations" );
818 0 : head.append<int>( "NITER", n, "actual number of iterations" );
819 0 : head.append<std::string>( "XREASON", reason, "reason for convergence" );
820 0 : head.append<realT>( "FPMTRANS", trans, "transmission of FPM" );
821 :
822 0 : fits::fitsFile<double> ff;
823 :
824 0 : std::string fname = "coron/" + cname + "_apod.fits";
825 :
826 0 : ff.write( fname, m_pupilApodizer, head );
827 :
828 0 : fname = "coron/" + cname + "_fpm.fits";
829 0 : improc::eigenImage<fpmaskFloatT> fpm( m_focalMask.rows(), m_focalMask.cols() );
830 0 : for( int r = 0; r < m_focalMask.rows(); ++r )
831 : {
832 0 : for( int c = 0; c < m_focalMask.cols(); ++c )
833 : {
834 0 : fpm( r, c ) = m_focalMask( r, c ).real();
835 : }
836 : }
837 :
838 0 : ff.write( fname, fpm, head );
839 :
840 0 : fname = "coron/" + cname + "_lyot.fits";
841 0 : ff.write( fname, m_lyotStop, head );
842 0 : }
843 :
844 : } // namespace wfp
845 : } // namespace mx
846 :
847 : #endif //__lyotCoronagraph_hpp__
|
