add pupil_vlt

maoppy/instrument.py

@@ -15,6 +15,7 @@ from scipy.interpolate import interp2d
 
 from maoppy.utils import circarr as _circarr
 from maoppy.utils import RAD2ARCSEC as _RAD2ARCSEC
+from maoppy.utils import binning as _binning
 
 #%% INSTRUMENT CLASS
 class Instrument:
@@ -131,6 +132,65 @@ class Instrument:
         return samp*(self.resolution_rad*self.D)
 
 
+def pupil_vlt(shape, angle=0, oversamp=1):
+    """
+    Returns the VLT pupil, rotated by <angle> (rad).
+    The keyword <oversamp> allows better sampling for smoother spiders
+    """
+    ### VLT PARAMETERS
+    occ = 0.14
+    D = 8
+    spider_width_m = 0.05
+    ### NUMERICAL VALUES
+    sx = shape[0] * oversamp
+    sy = shape[1] * oversamp
+    Rpix = min([sx,sy])/2 # same convention as 'Instrument.pupil'
+    pix_size_m = D/(2*Rpix)
+    spider_width_pix = spider_width_m / pix_size_m
+    ### CENTERED and ROTATED ARRAYS
+    xx0,yy0 = np.mgrid[0:sx,0:sy]
+    cx = (sx-1) / 2 # same convention as 'circarr'
+    cy = (sy-1) / 2
+    xx0 = xx0 - cx
+    yy0 = yy0 - cy
+    xx = np.cos(angle)*xx0 + np.sin(angle)*yy0
+    yy = np.cos(angle)*yy0 - np.sin(angle)*xx0
+    ### CIRCULAR PUPIL (with central occultation)
+    rr = np.sqrt(xx**2 + yy**2)
+    pup = (rr < Rpix) * (rr >= (Rpix*occ))
+    ### SPIDERS 1 and 2
+    x0 = - Rpix
+    y0 = 0
+    x1 = - Rpix*occ/np.sqrt(2)
+    y1 = Rpix*occ/np.sqrt(2)
+    a = (y1-y0)/(x1-x0)
+    spider1 = np.abs(yy-(a*(xx-x1)+y1)) <= (spider_width_pix/2)
+    spider1 *= (xx<=x1)
+    x0 = 0
+    y0 = - Rpix
+    a = -(x1-x0)/(y1-y0)
+    spider2 = np.abs(xx-(a*(yy-y1)+x1)) <= (spider_width_pix/2)
+    spider2 *= (yy>=y1)
+    ### SPIDERS 3 and 4
+    x0 = Rpix
+    y0 = 0
+    x1 = Rpix*occ/np.sqrt(2)
+    y1 = - Rpix*occ/np.sqrt(2)
+    a = (y1-y0)/(x1-x0)
+    spider3 = np.abs(yy-(a*(xx-x1)+y1)) <= (spider_width_pix/2)
+    spider3 *= (xx>=x1)
+    x0 = 0
+    y0 = Rpix
+    a = -(x1-x0)/(y1-y0)
+    spider4 = np.abs(xx-(a*(yy-y1)+x1)) <= (spider_width_pix/2)
+    spider4 *= (yy<=y1)
+    ### FINALIZE
+    pup_total = pup*(1-spider1)*(1-spider2)*(1-spider3)*(1-spider4)
+    if oversamp>1:
+        return _binning(pup_total, oversamp) / oversamp**2
+    return pup_total
+
+
 #%% LOAD INSTRUMENT INSTANCES (make them attributes of this module)
 # this might be clumsy, should I make something like this:
 # from maoppy.instrument import load_instrument