39 import astropy.io.fits
as pfits
52 """ Initialize the coronagraph
54 wavelength_0 = p_corono._wavelength_0
55 delta_wav = p_corono._delta_wav
56 nb_wav = p_corono._nb_wav
60 p_corono.set_nb_wav(1)
61 p_corono.set_wav_vec(np.array([wavelength_0]))
63 p_corono.set_delta_wav(0.)
64 p_corono.set_wav_vec(np.array([wavelength_0]))
66 wav_vec = np.linspace(wavelength_0 - delta_wav / 2,
67 wavelength_0 + delta_wav / 2,
70 p_corono.set_wav_vec(wav_vec)
73 if p_corono._type == scons.CoronoType.SPHERE_APLC:
75 elif p_corono._type == scons.CoronoType.PERFECT:
83 """ Dedicated function for SPHERE APLC coronagraph init
86 APLC data : SPHERE user manual, appendix A6 'NIR coronagraphs'
87 https://www.eso.org/sci/facilities/paranal/instruments/sphere/doc.html
88 IRDIS data : ESO instrument description
89 https://www.eso.org/sci/facilities/paranal/instruments/sphere/inst.html
92 p_corono.set_apodizer_name(scons.ApodizerType.SPHERE_APLC_APO1)
96 if p_corono._focal_plane_mask_name ==
None:
97 p_corono.set_focal_plane_mask_name(scons.FpmType.SPHERE_APLC_fpm_ALC2)
101 p_corono.set_lyot_stop_name(scons.LyotStopType.SPHERE_APLC_LYOT_STOP)
105 if p_corono._dim_image ==
None:
106 p_corono.set_dim_image(256)
107 if p_corono._image_sampling ==
None:
108 irdis_plate_scale = 12.25
109 VLT_pupil_diameter = 8
110 lambda_over_D = p_corono._wavelength_0 * 1e-6 / VLT_pupil_diameter
111 image_sampling = (lambda_over_D * 180 / np.pi * 3600 * 1000) / irdis_plate_scale
112 p_corono.set_image_sampling(image_sampling)
115 """ Dedicated function for perfect coronagraph init
122 if p_corono._apodizer_name == scons.ApodizerType.SPHERE_APLC_APO1:
123 apodizer = util.make_sphere_apodizer(pupdiam)
124 elif p_corono._apodizer_name ==
None:
125 apodizer = np.ones((pupdiam, pupdiam))
126 elif isinstance(p_corono._apodizer_name, str):
127 if not os.path.exists(p_corono._apodizer_name):
128 error_string =
"apodizer keyword (or path) '{}'".format(p_corono._apodizer_name) \
129 +
" is not a known keyword (or path)"
130 raise ValueError(error_string)
131 apodizer = pfits.getdata(p_corono._apodizer_name)
133 raise TypeError(
'apodizer name should be a string')
134 p_corono.set_apodizer(apodizer)
137 """ Focal plane mask init
139 if p_corono._focal_plane_mask_name == scons.FpmType.CLASSICAL_LYOT:
140 classical_lyot =
True
141 elif p_corono._focal_plane_mask_name
in (scons.FpmType.SPHERE_APLC_fpm_ALC1,
142 scons.FpmType.SPHERE_APLC_fpm_ALC2,
143 scons.FpmType.SPHERE_APLC_fpm_ALC3):
144 classical_lyot =
True
145 if (p_corono._focal_plane_mask_name == scons.FpmType.SPHERE_APLC_fpm_ALC1):
146 fpm_radius_in_mas = 145 / 2
147 elif (p_corono._focal_plane_mask_name == scons.FpmType.SPHERE_APLC_fpm_ALC2):
148 fpm_radius_in_mas = 185 / 2
149 elif (p_corono._focal_plane_mask_name == scons.FpmType.SPHERE_APLC_fpm_ALC3):
150 fpm_radius_in_mas = 240 / 2
151 VLT_pupil_diameter = 8
152 lambda_over_D = p_corono._wavelength_0 * 1e-6 / VLT_pupil_diameter
153 fpm_radius = fpm_radius_in_mas / (lambda_over_D * 180 / np.pi * 3600 * 1000)
154 p_corono.set_lyot_fpm_radius(fpm_radius)
156 classical_lyot =
False
159 p_corono.set_babinet_trick(
True)
160 if p_corono._fpm_sampling ==
None:
161 p_corono.set_fpm_sampling(20.)
162 lyot_fpm_radius_in_pix = p_corono._fpm_sampling * p_corono._lyot_fpm_radius
163 dim_fpm = 2 * int(lyot_fpm_radius_in_pix) + 2
164 p_corono.set_dim_fpm(dim_fpm)
165 fpm = util.classical_lyot_fpm(p_corono._lyot_fpm_radius,
167 p_corono._fpm_sampling,
170 elif isinstance(p_corono._focal_plane_mask_name, str):
171 if not os.path.exists(p_corono._focal_plane_mask_name):
172 error_string =
"focal plane mask keyword (or path) '{}'".format(p_corono._focal_plane_mask_name) \
173 +
" is not a known keyword (or path)"
174 raise ValueError(error_string)
175 fpm_array = pfits.getdata(p_corono._focal_plane_mask_name)
176 p_corono.set_dim_fpm(fpm_array.shape[0])
177 print(p_corono._dim_fpm)
178 if p_corono._fpm_sampling ==
None:
179 p_corono.set_fpm_sampling(p_corono._image_sampling)
180 if len(fpm_array.shape) == 2:
181 fpm = [fpm_array] * p_corono._nb_wav
182 elif len(fpm_array.shape) == 3:
184 for i
in range(p_corono._nb_wav):
185 fpm.append(fpm_array[:, :, i])
186 p_corono.set_focal_plane_mask(fpm)
191 if p_corono._lyot_stop_name == scons.LyotStopType.SPHERE_APLC_LYOT_STOP:
192 lyot_stop = util.make_sphere_lyot_stop(pupdiam)
193 elif p_corono._lyot_stop_name ==
None:
194 lyot_stop = np.ones((pupdiam, pupdiam))
195 elif isinstance(p_corono._lyot_stop_name, str):
196 if not os.path.exists(p_corono._lyot_stop_name):
197 error_string =
"Lyot stop keyword (or path) '{}'".format(p_corono._lyot_stop_name) \
198 +
" is not a known keyword (or path)"
199 raise ValueError(error_string)
200 lyot_stop = pfits.getdata(p_corono._lyot_stop_name)
202 raise TypeError(
'Lyot stop name should be a string')
203 p_corono.set_lyot_stop(lyot_stop)
205 def init_mft(p_corono: conf.Param_corono, pupdiam, planes, center_on_pixel=
False):
206 """ Initialize mft matrices
208 dim_fpm = p_corono._dim_fpm
209 fpm_sampling = p_corono._fpm_sampling
210 dim_image = p_corono._dim_image
211 image_sampling = p_corono._image_sampling
212 wavelength_0 = p_corono._wavelength_0
213 wav_vec = p_corono._wav_vec
215 norm0 = np.zeros(len(wav_vec))
217 if planes ==
'apod_to_fpm':
218 AA_apod_to_fpm = np.zeros((dim_fpm, pupdiam, len(wav_vec)), dtype=np.complex64)
219 BB_apod_to_fpm = np.zeros((pupdiam, dim_fpm, len(wav_vec)), dtype=np.complex64)
221 for w, wavelength
in enumerate(wav_vec):
222 wav_ratio = wavelength / wavelength_0
223 nbres = dim_fpm / fpm_sampling / wav_ratio
224 AA_apod_to_fpm[:,:,w], BB_apod_to_fpm[:,:,w], norm0[w] =
mft_matrices(pupdiam,
227 return AA_apod_to_fpm, BB_apod_to_fpm, norm0
229 elif planes ==
'fpm_to_lyot':
230 AA_fpm_to_lyot = np.zeros((pupdiam, dim_fpm, len(wav_vec)), dtype=np.complex64)
231 BB_fpm_to_lyot = np.zeros((dim_fpm, pupdiam, len(wav_vec)), dtype=np.complex64)
233 for w, wavelength
in enumerate(wav_vec):
234 wav_ratio = wavelength / wavelength_0
235 nbres = dim_fpm / fpm_sampling / wav_ratio
237 AA_fpm_to_lyot[:,:,w], BB_fpm_to_lyot[:,:,w], norm0[w] =
mft_matrices(dim_fpm,
238 pupdiam, nbres,inverse=
True)
240 return AA_fpm_to_lyot, BB_fpm_to_lyot, norm0
242 elif planes ==
'lyot_to_image':
243 AA_lyot_to_image = np.zeros((dim_image, pupdiam, len(wav_vec)), dtype=np.complex64)
244 BB_lyot_to_image = np.zeros((pupdiam, dim_image, len(wav_vec)), dtype=np.complex64)
246 for w, wavelength
in enumerate(wav_vec):
247 wav_ratio = wavelength / wavelength_0
248 nbres = dim_image / image_sampling / wav_ratio
250 AA_lyot_to_image[:,:,w], BB_lyot_to_image[:,:,w], norm0[w] =
mft_matrices(pupdiam,
256 AA_lyot_to_image[:,:,w], BB_lyot_to_image[:,:,w], norm0[w] =
mft_matrices(pupdiam,
260 return AA_lyot_to_image, BB_lyot_to_image, norm0
282 error_string_dim_input =
"'dim_input' must be an int (square input) or tuple of ints of dimension 2"
283 if np.isscalar(dim_input):
284 if isinstance(dim_input, int):
285 dim_input_x = dim_input
286 dim_input_y = dim_input
288 raise TypeError(dim_input)
289 elif isinstance(dim_input, tuple):
290 if all(isinstance(dims, int)
for dims
in dim_input) & (len(dim_input) == 2):
291 dim_input_x = dim_input[0]
292 dim_input_y = dim_input[1]
294 raise TypeError(error_string_dim_input)
296 raise TypeError(error_string_dim_input)
299 if real_dim_input ==
None:
300 real_dim_input = dim_input
301 error_string_real_dim_input =
"'real_dim_input' must be an int (square input pupil) or tuple of ints of dimension 2"
302 if np.isscalar(real_dim_input):
303 if isinstance(real_dim_input, int):
304 real_dim_input_x = real_dim_input
305 real_dim_input_y = real_dim_input
307 raise TypeError(error_string_real_dim_input)
308 elif isinstance(real_dim_input, tuple):
309 if all(isinstance(dims, int)
for dims
in real_dim_input) & (len(real_dim_input) == 2):
310 real_dim_input_x = real_dim_input[0]
311 real_dim_input_y = real_dim_input[1]
313 raise TypeError(error_string_real_dim_input)
315 raise TypeError(error_string_real_dim_input)
318 error_string_dim_output =
"'dim_output' must be an int (square output) or tuple of ints of dimension 2"
319 if np.isscalar(dim_output):
320 if isinstance(dim_output, int):
321 dim_output_x = dim_output
322 dim_output_y = dim_output
324 raise TypeError(error_string_dim_output)
325 elif isinstance(dim_output, tuple):
326 if all(isinstance(dims, int)
for dims
in dim_output) & (len(dim_output) == 2):
327 dim_output_x = dim_output[0]
328 dim_output_y = dim_output[1]
330 raise TypeError(error_string_dim_output)
332 raise TypeError(error_string_dim_output)
335 error_string_nbr =
"'nbres' must be an float or int (square output) or tuple of float or int of dimension 2"
336 if np.isscalar(nbres):
337 if isinstance(nbres, (float, int)):
338 nbresx = float(nbres)
339 nbresy = float(nbres)
341 raise TypeError(error_string_nbr)
342 elif isinstance(nbres, tuple):
343 if all(isinstance(nbresi, (float, int))
for nbresi
in nbres) & (len(nbres) == 2):
344 nbresx = float(nbres[0])
345 nbresy = float(nbres[1])
347 raise TypeError(error_string_nbr)
349 raise TypeError(error_string_nbr)
351 if real_dim_input
is not None:
352 nbresx = nbresx * dim_input_x / real_dim_input_x
353 nbresy = nbresy * dim_input_y / real_dim_input_y
355 X0 = dim_input_x / 2 + X_offset_input
356 Y0 = dim_input_y / 2 + Y_offset_input
358 X1 = dim_output_x / 2 + X_offset_output
359 Y1 = dim_output_y / 2 + Y_offset_output
361 xx0 = ((np.arange(dim_input_x) - X0 + 1 / 2) / dim_input_x)
362 xx1 = ((np.arange(dim_input_y) - Y0 + 1 / 2) / dim_input_y)
363 uu0 = ((np.arange(dim_output_x) - X1 + 1 / 2) / dim_output_x) * nbresx
364 uu1 = ((np.arange(dim_output_y) - Y1 + 1 / 2) / dim_output_y) * nbresy
367 if norm ==
'backward':
369 elif norm ==
'forward':
370 norm0 = nbresx * nbresy / dim_input_x / dim_input_y / dim_output_x / dim_output_y
371 elif norm ==
'ortho':
372 norm0 = np.sqrt(nbresx * nbresy / dim_input_x / dim_input_y / dim_output_x / dim_output_y)
373 sign_exponential = -1
376 if norm ==
'backward':
377 norm0 = nbresx * nbresy / dim_input_x / dim_input_y / dim_output_x / dim_output_y
378 elif norm ==
'forward':
380 elif norm ==
'ortho':
381 norm0 = np.sqrt(nbresx * nbresy / dim_input_x / dim_input_y / dim_output_x / dim_output_y)
384 AA = np.exp(sign_exponential * 1j * 2 * np.pi * np.outer(uu0, xx0)).astype(
'complex64')
385 BB = np.exp(sign_exponential * 1j * 2 * np.pi * np.outer(xx1, uu1)).astype(
'complex64')
390 """ Computes matrix multiplication for MFT
392 return norm * ((AA @ image.astype(
'complex64')) @ BB)
Parameter classes for COMPASS.
Numerical constants for shesha and config enumerations for safe-typing.
def init_apodizer(conf.Param_corono p_corono, pupdiam)
Apodizer init.
def init_focal_plane_mask(conf.Param_corono p_corono)
Focal plane mask init.
def init_perfect_coronagraph(conf.Param_corono p_corono, pupdiam)
Dedicated function for perfect coronagraph init.
def init_lyot_stop(conf.Param_corono p_corono, pupdiam)
Lyot stop init.
def mft_matrices(dim_input, dim_output, nbres, real_dim_input=None, inverse=False, norm='ortho', X_offset_input=0, Y_offset_input=0, X_offset_output=0, Y_offset_output=0)
docstring
def init_mft(conf.Param_corono p_corono, pupdiam, planes, center_on_pixel=False)
Initialize mft matrices.
def init_sphere_aplc(conf.Param_corono p_corono, pupdiam)
Dedicated function for SPHERE APLC coronagraph init.
def mft_multiplication(image, AA, BB, norm)
Computes matrix multiplication for MFT.
def init_coronagraph(conf.Param_corono p_corono, pupdiam)
Initialize the coronagraph.