package organization

README.md

@@ -29,10 +29,15 @@ Laboratoire d'Astrophysique de Marseille
 
 ## Reference
 
-If you may use the `Psfao` model of this library, the reference is
+The `Psfao` model of this library is fully described in the publication:
 
 [Fétick et al., August 2019, Astronomy and Astrophysics, Vol.628](https://www.aanda.org/articles/aa/abs/2019/08/aa35830-19/aa35830-19.html)
 
+For further reading on the `Psfao` model:
+
+[Beltramo-Martin et al., Astronomy and Astrophysics (forthcoming)](https://www.aanda.org/component/article?access=doi&doi=10.1051/0004-6361/202038679)
+
+
 ## License
 
 See the [LICENSE file](LICENSE.md)

maoppy/utils.py

@@ -6,13 +6,12 @@ Created on Mon May 27 17:27:44 2019
 @author: rfetick
 """
 
-import numpy as np
-from scipy.special import jv
-from numpy.fft import fft2, fftshift
+import numpy as _np
+from scipy.special import jv as _jv
 
 #%% CONSTANTS
 # Radian to arcsec
-RAD2ARCSEC = 180./np.pi * 3600.
+RAD2ARCSEC = 180./_np.pi * 3600.
 
 #%% IMAGE PROCESSING FUNCTIONS
 def circarr(Npix, center=(None, None)):
@@ -41,24 +40,24 @@ def circarr(Npix, center=(None, None)):
     if len(center) != 2: raise ValueError("Keyword `center` should be a tuple of 2 elements")
     if center[0] is None: center[0] = (Npix[0]-1) / 2
     if center[1] is None: center[1] = (Npix[1]-1) / 2
-    x, y = np.ogrid[0:Npix[0], 0:Npix[1]]
+    x, y = _np.ogrid[0:Npix[0], 0:Npix[1]]
     
-    return np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2)
+    return _np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2)
 
 def binning(image, k):
     """Bin an image by a factor `k`
     
     Example
     -------
-    >>> x,y = np.mgrid[0:10,0:10]
+    >>> x,y = _np.mgrid[0:10,0:10]
     >>> data = (-1)**x * (-1)**y
     >>> data_bin = binning(data,2)
     
     """
-    S = np.shape(image)
+    S = _np.shape(image)
     S0 = int(S[0] / k)
     S1 = int(S[1] / k)
-    out = np.zeros((S0, S1))
+    out = _np.zeros((S0, S1))
     for i in range(k):
         for j in range(k):
             out += image[i:k*S0:k, j:k*S1:k]
@@ -90,12 +89,12 @@ def airy(Npix, samp, occ=0., center=(None,None)):
     if samp <= 0: raise ValueError("You should ensure samp > 0")
     if occ < 0 or occ >= 1: raise ValueError("You should ensure 0 <= occ < 1")
 
-    r = circarr(Npix,center=center) * np.pi / samp
-    index = np.where(r == 0)
+    r = circarr(Npix,center=center) * _np.pi / samp
+    index = _np.where(r == 0)
     r[index] = 1.0
     
-    if occ > 0: a = 4.0/((1-occ**2)**2) * ((jv(1,r) - occ*jv(1,occ*r))/r)**2
-    else: a = 4.0 * (jv(1,r)/r) ** 2
+    if occ > 0: a = 4.0/((1-occ**2)**2) * ((_jv(1,r) - occ*_jv(1,occ*r))/r)**2
+    else: a = 4.0 * (_jv(1,r)/r) ** 2
 
     a[index] = 1.0
     return a / a.sum()
@@ -130,7 +129,7 @@ def imcenter(im, size, maxi=True, GC=False, center=None):
     ### SET to False the default choice, since user wants another option
     if (GC == True) or center != None: maxi = False
 
-    sX,sY = np.shape(im)
+    sX,sY = _np.shape(im)
 
     if len(size) != 2: raise ValueError("size keyword should contain 2 elements")
     if (size[0] < 0) or (size[1] < 0): raise ValueError("size keyword should contain only positive numbers")
@@ -140,23 +139,23 @@ def imcenter(im, size, maxi=True, GC=False, center=None):
 
     ### CENTER ON MAX, default method
     if maxi:
-        index = np.where(im == im.max())
+        index = _np.where(im == im.max())
         if len(index[0]) > 1: print("Imcenter warning: more than one maximum found")
         cX = index[0][0]
         cY = index[1][0]
 
     ### CENTER ON GRAVITY CENTER
     if GC:
-        x = np.arange(sX)
-        y = np.arange(sY)
-        cX = int(np.round(((x * im.sum(axis=1)).sum()) / (im.sum())))
-        cY = int(np.round(((y * im.sum(axis=0)).sum()) / (im.sum())))
+        x = _np.arange(sX)
+        y = _np.arange(sY)
+        cX = int(_np.round(((x * im.sum(axis=1)).sum()) / (im.sum())))
+        cY = int(_np.round(((y * im.sum(axis=0)).sum()) / (im.sum())))
 
     ### COMMON BLOCK FOR CENTERING!!!
     if (cX - size[0] / 2 < 0) or (cY - size[1] / 2 < 0) or (cX + size[0] / 2 >= sX) or (cY + size[1] / 2 >= sY):
         raise ValueError("Defined size is too big, cannot center image")
-    vX = np.arange(cX - size[0] / 2, cX + size[0] / 2, dtype=int)
-    vY = np.arange(cY - size[1] / 2, cY + size[1] / 2, dtype=int)
+    vX = _np.arange(cX - size[0] / 2, cX + size[0] / 2, dtype=int)
+    vY = _np.arange(cY - size[1] / 2, cY + size[1] / 2, dtype=int)
     newim = im[vX][:, vY]
 
     return newim
@@ -183,9 +182,9 @@ def circavg(tab, center=(None, None)):
     if tab.ndim != 2: raise ValueError("Input `tab` should be a 2D array")
     nx, ny = tab.shape
     r = circarr((int(nx), int(ny)), center=center)
-    cAVG = np.zeros(int(r.max()))
+    cAVG = _np.zeros(int(r.max()))
     for i in range(int(r.max())):
-        index = np.where((r >= i) * (r < (i + 1)))
+        index = _np.where((r >= i) * (r < (i + 1)))
         cAVG[i] = tab[index[0], index[1]].sum() / index[0].size
     return cAVG
 
@@ -196,7 +195,7 @@ def circavgplt(tab, center=(None, None), dtype=float):
         `y` contains the circular average of the 2D `tab`
     Output is symmetrized with respect to x=0 for a symmetric plot"""
     y = circavg(tab, center=center)
-    x = np.arange(len(y), dtype=dtype)
+    x = _np.arange(len(y), dtype=dtype)
     x2 = -x[:0:-1]
     y2 = y[:0:-1]
-    return (np.concatenate((x2, x)), np.concatenate((y2, y)))
+    return (_np.concatenate((x2, x)), _np.concatenate((y2, y)))

setup.py

@@ -18,5 +18,5 @@ setup(name='MAOPPY',
       package_data={'maoppy': ['data/*.yml']},
       include_package_data=True,
       packages=find_packages(exclude=['example','test']),
-      requires=['numpy','scipy','astropy'],
+      requires=['numpy','scipy','astropy','pyyaml'],
       zip_safe=False)