Merge branch 'develop' into master

README.md

@@ -40,8 +40,9 @@ If you use this software, or part it, please reference the code from the Astroph
 
 - [Vigan & N'Diaye, 2018, ascl:1806.003](https://ui.adsabs.harvard.edu/abs/2018ascl.soft06003V/abstract)
 
-Authors
--------
+Authors and contributors
+------------------------
 
 - Arthur Vigan, [[email protected]](mailto:[email protected])
 - Mamadou N'Diaye, [[email protected]](mailto:[email protected])
+- Raphaël Pourcelot, [[email protected]](mailto:[email protected])

examples/advanced_pyzelda_analysis.py

@@ -0,0 +1,101 @@
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+import numpy as np
+
+import pyzelda.zelda as zelda
+import pyzelda.ztools as ztools
+import pyzelda.utils.aperture as aperture
+import pyzelda.utils.zernike as zernike
+import pyzelda.utils.mft as mft
+
+from pathlib import Path
+
+#%% parameters
+
+# path = Path('/Users/mndiaye/Dropbox/python/zelda/pyZELDA/')
+path = Path('/Users/avigan/Work/GitHub/pyZELDA/data/')
+# path = Path('D:/Programmes/GitHub/pyZELDA/')
+
+wave = 1.642e-6
+
+z = zelda.Sensor('CLASSIC')
+
+#%%################################
+# Creation of apodized pupil data #
+###################################
+
+# Aberration: astigmatism
+basis = np.nan_to_num(zernike.zernike_basis(nterms=5, npix=z.pupil_diameter))*1e-9
+aberration = 20*basis[4]
+gaussian_std = 1
+
+# Gaussian apodized entrance pupil
+x_vals = np.arange(z.pupil_diameter)
+xx, yy = np.meshgrid(x_vals, x_vals)
+cx, cy = x_vals[z.pupil_diameter//2], x_vals[z.pupil_diameter//2]
+r = np.sqrt((xx-cx)**2+(yy-cy)**2)
+r = r/r[0, z.pupil_diameter//2]
+
+apodizer = np.exp(-(r**2)/2/gaussian_std**2)
+
+#%%#############################
+# Simulation of the wavefront. #
+################################
+
+# Clear pupil image
+clear_pupil = abs(apodizer)**2
+
+zelda_pupil = ztools.propagate_opd_map(aberration, z.mask_diameter, z.mask_depth, z.mask_substrate,
+                                        z.mask_Fratio, z.pupil_diameter, apodizer, wave)
+
+#%%#########################
+# Wavefront reconstruction #
+############################
+
+# Standard analysis
+z_opd_standard = z.analyze(clear_pupil, zelda_pupil, wave)
+
+# Advanced analysis
+pupil_roi = aperture.disc(z.pupil_diameter, z.pupil_diameter, diameter=True, cpix=False)
+z_opd_advanced = z.analyze(clear_pupil, zelda_pupil, wave,
+                           use_arbitrary_amplitude=True,
+                           refwave_from_clear=True,
+                           cpix=False, pupil_roi=pupil_roi)
+
+#%% plot
+
+fig = plt.figure(0, figsize=(24, 4))
+plt.clf()
+
+gs = gridspec.GridSpec(ncols=7, nrows=1, figure=fig, width_ratios=[.1,1,1,1,1,1,.1])
+
+ax = fig.add_subplot(gs[0,1])
+mappable = ax.imshow(clear_pupil, aspect='equal', vmin=0, vmax=1)
+ax.set_title('Clear pupil')
+
+ax = fig.add_subplot(gs[0,0])
+cbar1 = fig.colorbar(mappable=mappable, cax=ax)
+cbar1.set_label('Normalized intensity')
+
+ax = fig.add_subplot(gs[0,2])
+ax.imshow(zelda_pupil, aspect='equal', vmin=0, vmax=1)
+ax.set_title('ZELDA pupil')
+
+ax = fig.add_subplot(gs[0,3])
+ax.imshow(aberration*1e9, aspect='equal', vmin=-40, vmax=40, cmap='magma')
+ax.set_title('Introduced aberration (nm)')
+
+ax = fig.add_subplot(gs[0,4])
+cax = ax.imshow(z_opd_standard, aspect='equal', vmin=-40, vmax=40, cmap='magma')
+ax.set_title('Reconstructed aberration - standard')
+
+ax = fig.add_subplot(gs[0,5])
+cax = ax.imshow(z_opd_advanced, aspect='equal', vmin=-40, vmax=40, cmap='magma')
+ax.set_title('Reconstructed aberration - advanced')
+
+ax = fig.add_subplot(gs[0,6])
+cbar = fig.colorbar(mappable=cax, cax=ax)
+cbar.set_label('OPD [nm]')
+
+plt.tight_layout()
+plt.show()

examples/simple_zelda_analysis.py

@@ -7,6 +7,8 @@
 
 from pathlib import Path
 
+#%% parameters
+
 # path = Path('/Users/mndiaye/Dropbox/python/zelda/pyZELDA/')
 path = Path('/Users/avigan/Work/GitHub/pyZELDA/data/')
 # path = Path('D:/Programmes/GitHub/pyZELDA/')
@@ -26,35 +28,37 @@
 dark_file  = 'SPHERE_GEN_IRDIS057_0003'
 pupil_tel  = True
 
-# ZELDA analysis
+#%% ZELDA analysis
 z = zelda.Sensor('SPHERE-IRDIS', pupil_telescope=pupil_tel)
 
 clear_pupil, zelda_pupil, center = z.read_files(path, clear_pupil_files, zelda_pupil_files, dark_file,
-                                                collapse_clear=True, collapse_zelda=True)
+                                                collapse_clear=True, collapse_zelda=False)
 
 opd_map = z.analyze(clear_pupil, zelda_pupil, wave=wave)
+if opd_map.ndim == 3:
+    opd_map = opd_map.mean(axis=0)
 
 # decomposition on Zernike polynomials
-basis, coeff, opd_zern = ztools.zernike_expand(opd_map.mean(axis=0), 20)
+basis, coeff, opd_zern = ztools.zernike_expand(opd_map, 20)
 
-# plot
+#%% plot
 fig = plt.figure(0, figsize=(16, 4))
 plt.clf()
 
 ax = fig.add_subplot(141)
-ax.imshow(clear_pupil.mean(axis=0), aspect='equal', vmin=0, vmax=15000, origin=1)
+ax.imshow(clear_pupil.mean(axis=0), aspect='equal', vmin=0, vmax=2000)
 ax.set_title('Clear pupil')
 
 ax = fig.add_subplot(142)
-ax.imshow(zelda_pupil.mean(axis=0), aspect='equal', vmin=0, vmax=15000, origin=1)
+ax.imshow(zelda_pupil.mean(axis=0), aspect='equal', vmin=0, vmax=2000)
 ax.set_title('ZELDA pupil')
 
 ax = fig.add_subplot(143)
-ax.imshow(opd_map.mean(axis=0), aspect='equal', vmin=-150, vmax=150, cmap='magma', origin=1)
+ax.imshow(opd_map, aspect='equal', vmin=-150, vmax=150, cmap='magma')
 ax.set_title('OPD map')
 
 ax = fig.add_subplot(144)
-cax = ax.imshow(opd_zern.mean(axis=0), aspect='equal', vmin=-150, vmax=150, cmap='magma', origin=1)
+cax = ax.imshow(opd_zern.mean(axis=0), aspect='equal', vmin=-150, vmax=150, cmap='magma')
 ax.set_title('Zernike projected OPD map')
 
 cbar = fig.colorbar(cax)

pyzelda/__init__.py

@@ -1,7 +1,7 @@
 __author__    = "Arthur Vigan, Mamadou N'Diaye"
 __copyright__ = "Copyright (C) 2017-2020 Arthur Vigan, Mamadou N'Diaye"
 __license__   = 'MIT'
-__version__   = '1.1'
+__version__   = '1.2'
 
 # direct access to sensor object
 from .zelda import Sensor

pyzelda/sphere/sequence.py

@@ -332,7 +332,7 @@ def read_info(root):
     return info_files, info_frames, info_frames_ref
 
 
-def process(root, sequence_type='temporal', correction_factor=1):
+def process(root, sequence_type='temporal', correction_factor=1, unit='m'):
     '''Process a complete sequence of ZELDA data
 
     The processing centers the data and performs the ZELDA analysis to
@@ -352,6 +352,8 @@ def process(root, sequence_type='temporal', correction_factor=1):
         Amplitude correction factor for the OPD maps. 
         Default is 1 (no correction)
 
+    unit : str
+        Unit for the processed cube. Can me either m, um or nm. Default is m
     '''
 
     # read info
@@ -367,7 +369,17 @@ def process(root, sequence_type='temporal', correction_factor=1):
 
     # create sensor
     z = zelda.Sensor('SPHERE-IRDIS')
-
+
+    # apply unit and correction factor
+    if unit == 'm':
+        pass
+    elif unit == 'um':
+        correction_factor *= 1e6
+    elif unit == 'nm':
+        correction_factor *= 1e9
+    else:
+        raise ValueError(f'Unknown output unit {unit}')
+
     # read and analyse
     print('ZELDA analysis')
     if sequence_type == 'temporal':
@@ -382,9 +394,8 @@ def process(root, sequence_type='temporal', correction_factor=1):
             # analyse
             opd_cube = z.analyze(clear_pupil, zelda_pupils, wave=1.642e-6)
 
-            # correction factor
-            if correction_factor != 1:
-                opd_cube *= correction_factor
+            # correction factor and unit
+            opd_cube *= correction_factor
 
             fits.writeto(os.path.join(root, 'processed', zelda_pupil_files[f]+'_opd.fits'), opd_cube, overwrite=True)
 
@@ -436,9 +447,8 @@ def process(root, sequence_type='temporal', correction_factor=1):
                 # analyse
                 opd_cube[img] = z.analyze(clear_pupil[img_ref], zelda_pupils[img], wave=1.642e-6)
 
-            # correction factor
-            if correction_factor != 1:
-                opd_cube *= correction_factor
+            # correction factor and unit
+            opd_cube *= correction_factor
 
             fits.writeto(os.path.join(root, 'processed', zelda_pupil_files[f]+'_opd.fits'), opd_cube, overwrite=True)
             del opd_cube