PSF is normalized to unit energy at infinity and not any more on the FoV

maoppy/psfmodel.py

@@ -473,6 +473,7 @@ class Psfao(ParametricPSF):
         Returns
         -------
         psd : numpy.array (dim=2)
+        integral : float : the integral of the `psd` up to infinity
             
         """
         self.check_parameters(x0)
@@ -500,7 +501,14 @@ class Psfao(ParametricPSF):
         PSD += (F2 < Fao**2.) * np.abs(C + A*moffat(f2D,param,norm=Fao))
         # Set PSD = 0 at null frequency
         PSD[Nx_over//2,Ny_over//2] = 0.0
-        return PSD
+        
+        # Compute PSD integral up to infinity
+        fmax = np.min([Nx_over//2,Ny_over//2])*pix2freq
+        integral_in = np.sum(PSD*(F2<(fmax**2))) * pix2freq**2 # numerical sum
+        integral_out = 0.0229*6*np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
+        integral = integral_in + integral_out
+        
+        return PSD, integral
     
     def otf(self,x0,dx=0,dy=0,_caller='user'):
         """
@@ -520,15 +528,15 @@ class Psfao(ParametricPSF):
         return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
     
     def _otf_turbulent(self,x0):
-        PSD = self.psd(x0)
+        PSD,integral = self.psd(x0)
         L = self.system.D * self._samp_over
         Bg = fft2(fftshift(PSD)) / L**2
-        Dphi = fftshift(np.real(2 * (Bg[0, 0] - Bg)))
+        #Dphi = fftshift(np.real(2 * (Bg[0, 0] - Bg))) # normalized on the numerical FoV 
+        Dphi = fftshift(np.real(2 * (integral - Bg))) # normalized up to infinity
         return np.exp(-Dphi/2.) 
     
     @lru_cache(maxsize=2)
     def _otf_diffraction(self):
-        #TODO: @lru_cache to prevent unecessary calls to this? (2 FFT)
         Nx_over = self.Npix[0]*self._k
         Ny_over = self.Npix[1]*self._k
         
@@ -566,7 +574,6 @@ class Psfao(ParametricPSF):
         """
         
         out = np.real(fftshift(ifft2(fftshift(self.otf(x0,dx=dx,dy=dy,_caller='self')))))
-        out = out/out.sum() # ensure unit energy on the field of view
         
         if self._k==1:
             return out

maoppy/test/test_psfmodel.py

@@ -82,21 +82,21 @@ class TestPsfao(unittest.TestCase):
         theta = 0
         beta = 1.5
         # Integral of the halo
-        psd = P.psd([r0,C,A,alpha,ellip,theta,beta])
+        psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
         int_num = np.sum(psd)*df*df
         int_ana = 0.023 * 6*np.pi/5 * (fao*r0)**(-5.0/3.0)
         self.assertAlmostEqual(int_num,int_ana,delta=1e-2)
         # Integral of the constant
         r0 = np.inf
         C = 1e-2
-        psd = P.psd([r0,C,A,alpha,ellip,theta,beta])
+        psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
         int_num = np.sum(psd)*df*df
         int_ana = C*np.pi*fao**2.0 - df*df # remove central pixel
         self.assertAlmostEqual(int_num,int_ana,delta=1e-2)
         # Integral of the Moffat
         C = 0.0
         A = 1.0
-        psd = P.psd([r0,C,A,alpha,ellip,theta,beta])
+        psd,_ = P.psd([r0,C,A,alpha,ellip,theta,beta])
         int_num = np.sum(psd)*df*df
         int_ana = A
         self.assertAlmostEqual(int_num,int_ana,delta=1e-2)