Remove central pixel for Moffat energy computation

maoppy/example/fit_muse_simulated_psf_small_fov.py

@@ -33,7 +33,7 @@ wvl = 600*1e-9 # wavelength [m]
 
 flux = 1e6
 ron = muse_nfm.ron
-bckgd = 1.0
+bckgd = 10.0
 
 #%% Initialize PSF model
 samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
@@ -56,7 +56,7 @@ noise = ron*np.random.randn(npix,npix)
 image = flux*psf + bckgd + noise
 
 #%% Crop the image
-ncrop = 64
+ncrop = 32
 im_crop = image[npix//2-ncrop//2:npix//2+ncrop//2,
                 npix//2-ncrop//2:npix//2+ncrop//2]
 

maoppy/psfmodel.py

@@ -203,7 +203,7 @@ def reduced_coord(XY,ax,ay,theta,cx,cy):
     u = Rxx*(XY[0]-cx)**2 + Rxy*(XY[0]-cx)*(XY[1]-cy) + Ryy*(XY[1]-cy)**2
     return u    
 
-def moffat(XY, param, norm=None):
+def moffat(XY, param, norm=None, removeInside=0):
     """
     Compute a Moffat function on a meshgrid
     moff = Enorm * (1+u)^(-beta)
@@ -231,6 +231,9 @@ def moffat(XY, param, norm=None):
                     Enorm = (beta-1)/(pi*ax*ay)
         float,int - Energy normalization up to the radius defined by this value
                     Enorm = (beta-1)/(pi*ax*ay)*(1-(1+(R**2)/(ax*ay))**(1-beta))    
+    
+    removeInside: float (default=0)
+        Used to remove the central pixel in energy computation
     """
     if len(param)!=6: raise ValueError("Parameter `param` must contain exactly 6 elements, but input has %u elements"%len(param))
     ax,ay,theta,beta,cx,cy = param
@@ -239,13 +242,13 @@ def moffat(XY, param, norm=None):
     
     if norm is None:
         Enorm = 1
-    elif norm == _np.inf:
-        if beta<=1: raise ValueError("Cannot compute Moffat energy for beta<=1")
-        Enorm = (beta-1) / (_np.pi*ax*ay)
-    else:
+    else: # norm can be float or np.inf
+        if (beta<=1) and (norm==_np.inf): raise ValueError("Cannot compute Moffat energy for beta<=1")
         if beta==1: raise ValueError("Energy computation for beta=1.0 not implemented yet. Sorry!")
         Enorm = (beta-1) / (_np.pi*ax*ay)
-        Enorm = Enorm / (1 - (1 + (norm**2) / (ax*ay))**(1-beta))
+        Fout = (1 - (1 + (norm**2) / (ax*ay))**(1-beta))
+        Fin = (1 - (1 + (removeInside**2) / (ax*ay))**(1-beta))
+        Enorm = Enorm / (Fout-Fin)
     
     return Enorm * (1.0+u)**(-beta)
 
@@ -463,11 +466,13 @@ class ParametricPSFfromPSD(ParametricPSF):
         return self._xyarray * self._pix2freq
     
     def strehlOTF(self,x0):
+        """Compute the Strehl ratio based on the sum of the OTF"""
         otf = self.otf(x0,_caller='self')
-        otf_diff = self._otf_diffraction()
+        otf_diff = self._otfDiffraction()
         return _np.real(_np.sum(otf)/_np.sum(otf_diff))
     
     def strehlMarechal(self,x0):
+        """Compute the Strehl ratio based on the Maréchal approximation"""
         _,sig2 = self.psd(x0)
         return _np.exp(-sig2)
     
@@ -494,13 +499,13 @@ class ParametricPSFfromPSD(ParametricPSF):
         if (self._k > 1) and (_caller != 'self'):
             raise ValueError("Cannot call `otf(...)` when undersampled (functionality not implemented yet)")
         
-        OTF_TURBULENT = self._otf_turbulent(x0)
-        OTF_DIFFRACTION = self._otf_diffraction()
-        OTF_SHIFT = self._otf_shift(dx,dy)
+        otfTurbulent = self._otfTurbulent(x0)
+        otfDiffraction = self._otfDiffraction()
+        otfShift = self._otfShift(dx,dy)
         
-        return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
+        return otfTurbulent * otfDiffraction * otfShift
     
-    def _otf_turbulent(self,x0):
+    def _otfTurbulent(self,x0):
         PSD,integral = self.psd(x0)
         L = self.system.D * self._samp_over
         Bg = _fft.fft2(_fft.fftshift(PSD)) / L**2
@@ -509,7 +514,7 @@ class ParametricPSFfromPSD(ParametricPSF):
         return _np.exp(-Dphi/2.)
     
     #@_lru_cache(maxsize=2)
-    def _otf_diffraction(self):
+    def _otfDiffraction(self):
         nx, ny = self._nxnyOver
         NpupX = _np.ceil(nx/self._samp_over)
         NpupY = _np.ceil(ny/self._samp_over)
@@ -518,7 +523,7 @@ class ParametricPSFfromPSD(ParametricPSF):
         return _fft.fftshift(_np.abs(_fft.ifft2(_np.abs(_fft.fft2(tab)) ** 2)) / _np.sum(tab))
     
     #@_lru_cache(maxsize=2)
-    def _otf_shift(self,dx,dy):
+    def _otfShift(self,dx,dy):
         Y, X = self._xyarray
         ny,nx = X.shape
         return _np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny))
@@ -603,8 +608,8 @@ class Turbulent(ParametricPSFfromPSD):
         
         # Compute PSD integral up to infinity
         fmax = _np.min([nx0,ny0])*self._pix2freq
-        integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum
-        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
+        integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum (in the array's tangent circle)
+        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum (outside the array's tangent circle)
         integral = integral_in + integral_out
         
         return PSD, integral
@@ -686,7 +691,7 @@ class Psfao(ParametricPSFfromPSD):
         ay = alpha/ratio
         param = (ax,ay,theta,beta,0,0)
         
-        moff = moffat(f2D,param,norm=Fao) * (F2 < Fao**2.)
+        moff = moffat(f2D,param,norm=Fao,removeInside=self._pix2freq/2) * (F2 < Fao**2.)
         moff[nx0,ny0] = 0.0 # Set Moffat PSD = 0 at null frequency
         moff = moff / (_np.sum(moff)*self._pix2freq**2)  # normalize moffat numerically to get correct A=sigma² in the AO zone
         # Warning: Moffat numerical normalization generates strehlOTF jump when "_k" is changed