Manage phasemask (not tested yet)

maoppy/instrument.py

@@ -9,9 +9,29 @@ Created on Mon May 27 17:30:51 2019
 from maoppy.utils import circarr, RAD2ARCSEC
 import yaml
 import os
+from astropy.io import fits
+import numpy as np
+from scipy.interpolate import interp2d
 
-#%% INSTRUMENT CLASS AND ITS SUBCLASSES
 
+def _get_data_folder():
+    folder = os.path.abspath(__file__)
+    folder = os.sep.join(folder.split(os.sep)[0:-1])+os.sep+'data'+os.sep
+    return folder
+
+def load(name):
+    """Load an instrument saved in the `data` folder of MAOPPY"""
+    folder = _get_data_folder()
+    with open(folder+name.lower()+".yml",'r') as f:
+        d = yaml.full_load(f)
+    instru = Instrument(D=d['d'],occ=d['occ'],res=d['res'],gain=d['gain'],ron=d['ron'],Nact=d['nact'])
+    instru.name = d['name']
+    instru.fullname = d['fullname']
+    instru.filters = d['filters']
+    instru.phasemask_path = d['phasemask_path']
+    return instru
+
+#%% INSTRUMENT CLASS
 class Instrument(object):
     """Represents an optical system (telescope to detector)
     
@@ -56,7 +76,9 @@ class Instrument(object):
         self.fullname = "MAOPPY Instrument" # human readable name
         
         # phasemask (not used yet)
+        self.phasemask_enable = False
         self.phasemask_path = None
+        self._phasemask = None
         
     def __repr__(self):
         s  = "---------------------------------\n" 
@@ -72,6 +94,7 @@ class Instrument(object):
         s += "\n"
         s += "Detector   : gain=%.1f e-/ADU\n"%self.gain
         s += "             RON =%.1f e-\n"%self.ron
+        s += "Phasemask  : %s\n"%self.phasemask_enable
         s += "---------------------------------\n" 
         return s
     
@@ -84,10 +107,24 @@ class Instrument(object):
         return self.resolution_rad * RAD2ARCSEC * 1e3
     
     def pupil(self,Npix,wvl=None,samp=None):
-        """Returns the 2D array of the pupil transmission function"""
+        """Returns the 2D array of the pupil transmission function (complex data)"""
         Dpix = min(Npix)/2
         pup = circarr(Npix)
-        return (pup < Dpix) * (pup >= Dpix*self.occ)
+        if self.phasemask_enable:
+            if self._phasemask is None:
+                if self.phasemask_path is None: raise ValueError('phasemask_path must be defined')
+                p = fits.open(self.phasemask_path)[0].data * 1e-9 # fits data in nm, converted here to meter
+                x = np.arange(p.shape[0])/p.shape[0]
+                y = np.arange(p.shape[1])/p.shape[1]
+                self._phasemask = interp2d(x,y,p)
+            x = np.arange(Npix[0])/Npix[0] - 0.5/Npix[0]
+            y = np.arange(Npix[1])/Npix[1] - 0.5/Npix[1]
+            wf = self._phasemask(x,y)
+            if wvl is None: wvl = self.wvl(samp) # samp must be defined if wvl is None
+            wf = np.exp(2j*np.pi/wvl*wf)
+        else:
+            wf = 1.0 + 0j # complex type for output array, even if real data
+        return (pup < Dpix) * (pup >= Dpix*self.occ) * wf
     
     def samp(self,wvl):
         """Returns sampling value for the given wavelength"""
@@ -113,21 +150,6 @@ class Instrument(object):
         with open(filename,'w') as f: yaml.dump(data,f)
         
 
-# TODO: manage phasemask
-
-def load(name):
-    """Load an instrument saved in the `data` folder of MAOPPY"""
-    folder = os.path.abspath(__file__)
-    folder = os.sep.join(folder.split(os.sep)[0:-1])+os.sep+'data'+os.sep
-    with open(folder+name.lower()+".yml",'r') as f:
-        d = yaml.full_load(f)
-    instru = Instrument(D=d['d'],occ=d['occ'],res=d['res'],gain=d['gain'],ron=d['ron'],Nact=d['nact'])
-    instru.name = d['name']
-    instru.fullname = d['fullname']
-    instru.filters = d['filters']
-    instru.phasemask_path = d['phasemask_path']
-    return instru
-
 #%% OLD WAY to LOAD INSTRUMENTS (now use the `load` function)
 
 #ZIMPOL = Instrument(D=8.,occ=0.14,res=30*1e-6/1768.,gain=10.5,ron=20.,Nact=40)