add PsfaoAliasing class

maoppy/data/muse_nfm.ini

@@ -13,6 +13,8 @@ occ = 0.14
 [ao]
 ; Number of effectively controlled modes across the diameter
 nact = 34
+; Aliasing factor between aliased variance and fitting variance
+aliasing = 0.4
 
 [camera]
 ; Resolution of one pixel [mas]

maoppy/data/muse_wfm.ini

@@ -13,6 +13,8 @@ occ = 0.14
 [ao]
 ; Number of effectively controlled modes across the diameter
 nact = 34
+; Aliasing factor between aliased variance and fitting variance
+aliasing = 0.4
 
 [camera]
 ; Resolution of one pixel [mas]

maoppy/data/zimpol.ini

@@ -13,6 +13,8 @@ occ = 0.14
 [ao]
 ; Number of effectively controlled modes across the diameter
 nact = 40
+; Aliasing factor between aliased variance and fitting variance
+aliasing = 0.4
 
 [camera]
 ; Resolution of one pixel [mas]

maoppy/instrument.py

@@ -55,6 +55,7 @@ class Instrument:
         self.occ = occ # occultation diameter ratio
         self.filters = {}
         self.Nact = Nact
+        self.aliasing = 0
         
         self._resolution_rad = res
         self.gain = gain
@@ -160,6 +161,7 @@ def load_ini(pth):
     occ = float(config['telescope']['occ'])
     # [ao]
     nact = int(config['ao']['nact'])
+    aliasing = float(config['ao']['aliasing'])
     # [camera]
     res_mas = float(config['camera']['res_mas'])
     # [filters]
@@ -175,6 +177,7 @@ def load_ini(pth):
     instru.name = tag
     instru.fullname = name
     instru.filters = filters
+    instru.aliasing = aliasing
     return instru
     
 

maoppy/psfmodel.py

@@ -12,6 +12,7 @@ from scipy.fft import fft2, ifft2, fftshift, ifftshift
 from astropy.io import fits
 from maoppy.utils import binning
 from maoppy.psfutils import moffat, gauss, oversample, moffat_center, make_fits_hdr
+import warnings
 
 __all__ = ["ParametricPSF", "ConstantPSF", "Moffat", "Gaussian",
            "ParametricPSFfromPSD", "Turbulent", "Psfao"]
@@ -643,6 +644,9 @@ class Psfao(ParametricPSFfromPSD):
         
         r0, bck, amp, alpha, ratio, theta, beta = parampsd
         
+        if alpha < pix:
+            warnings.warn("`alpha < df`. Everything works fine but the value of `amp` might not correspond exactly to the PSD integral.")
+        
         # Von-Karman
         psd = (r0**(-5./3.)) * self._vk
         
@@ -735,3 +739,45 @@ class Psfao(ParametricPSFfromPSD):
         
         hdu = fits.PrimaryHDU(psf, hdr)
         hdu.writeto(filename)
+
+
+
+class PsfaoAliasing(Psfao):
+    """Ignore background and use aliasing instead"""
+    
+    def psd(self, parampsd, *args, grad=False, fovnorm=False, **kwargs):
+        # parampsd can be called with a dictionary
+        if isinstance(parampsd,dict):
+            parampsd = self.dict2list(parampsd)
+        
+        r0 = parampsd[0]
+        bck = self.system.aliasing*0.0229*6/5* r0**(-5/3) * self.cutoffAOfreq**(-11/3)
+        dbck_dr0 = -5/3 * bck/r0
+        parampsd[1] = bck
+        
+        out = super().psd(parampsd, *args, grad=grad, **kwargs)
+        
+        if grad:
+            # I need to update the gradient and its integral!
+            psd, integral, gg, _ = out
+            gg[0,...] += self._maskin * dbck_dr0
+            gg[1,...] = 0
+            
+            integral_g = np.zeros(len(parampsd))
+            pix = self.pix2freq
+            
+            if fovnorm:
+                for i in range(len(parampsd)):
+                    integral_g[i] = np.sum(gg[i,...]) * pix**2 # numerical sum
+            else:
+                nx0,ny0 = self._nullFreqIndex
+                fmax = np.min([nx0,ny0])*pix
+                integral_out = 0.0229*6*np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum outside
+                for i in range(len(parampsd)):
+                    integral_g[i] = np.sum(gg[i,...]*(self._f2<(fmax**2))) * pix**2 # numerical sum
+                integral_g[0] += integral_out*(-5./3)/r0
+                
+            psd, integral, gg, integral_g
+                
+        return out
+        
\ No newline at end of file