precompute Psfao _xyarray

maoppy/psfmodel.py

@@ -402,11 +402,20 @@ class ParametricPSFfromPSD(ParametricPSF):
         if (npix[0]%2) or (npix[1]%2): raise ValueError("Each npix component must be even")
         if system is None: raise ValueError("Keyword `system` must be defined")
         if samp is None: raise ValueError("Keyword `samp` must be defined")
-        self.npix = npix
+        self._npix = npix # "_" to bypass the _xyarray update, that will be made with samp setter
+        self.samp = samp # also init _xyarray
         self.system = system
-        self.samp = samp
         self._nparam = nparam
         
+    @property
+    def npix(self):
+        return self._npix
+    
+    @npix.setter
+    def npix(self,value):
+        self._npix = value
+        self._computeXYarray()
+        
     @property
     def samp(self):
         return self._samp_over/self._k
@@ -415,32 +424,13 @@ class ParametricPSFfromPSD(ParametricPSF):
     def samp(self,value):
         # Manage cases of undersampling
         self._samp_over, self._k = oversample(value)
-
-    def check_parameters(self,x0):
-        bd, bu = self.bounds
-        x0 = _np.array(x0)
-        bd = _np.array(bd)
-        bu = _np.array(bu)
-        if len(x0)!=self._nparam: raise ValueError('len(x0) is different from length of bounds')
-        if _np.any(x0<bd): raise ValueError('Lower bounds are not respected')
-        if _np.any(x0>bu): raise ValueError('Upper bounds are not respected')
+        self._computeXYarray()
         
-    def otf(self,x0,dx=0,dy=0,_caller='user'):
-        """
-        See __call__ for input arguments
-        Warning: result of otf will be unconsistent if undersampled!!!
-        This issue is solved with oversampling + binning in __call__ but not here
-        For the moment, the `_caller` keyword prevents user to misuse otf
-        """
-        
-        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)
-        
-        return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
+    def _computeXYarray(self):
+        nx,ny = self._nxnyOver
+        self._xyarray = _np.mgrid[0:nx, 0:ny].astype(float)
+        self._xyarray[0] -= nx//2
+        self._xyarray[1] -= ny//2
     
     #@_lru_cache(maxsize=2)
     @property
@@ -470,12 +460,7 @@ class ParametricPSFfromPSD(ParametricPSF):
         """
         Return the array of frequencies (correctly sampled above Shannon-Nyquist)
         """
-        Nx_over, Ny_over = self._nxnyOver
-        f2D = _np.mgrid[0:Nx_over, 0:Ny_over].astype(float)
-        f2D[0] -= Nx_over//2
-        f2D[1] -= Ny_over//2
-        f2D *= self._pix2freq
-        return f2D
+        return self._xyarray * self._pix2freq
     
     def strehlOTF(self,x0):
         otf = self.otf(x0,_caller='self')
@@ -489,6 +474,32 @@ class ParametricPSFfromPSD(ParametricPSF):
     def psd(self,x0):
         raise ValueError("ParametricPSFfromPSD is not to be instantiated. the `psd` method must be override in the subclasses")
     
+    def check_parameters(self,x0):
+        bd, bu = self.bounds
+        x0 = _np.array(x0)
+        bd = _np.array(bd)
+        bu = _np.array(bu)
+        if len(x0)!=self._nparam: raise ValueError('len(x0) is different from length of bounds')
+        if _np.any(x0<bd): raise ValueError('Lower bounds are not respected')
+        if _np.any(x0>bu): raise ValueError('Upper bounds are not respected')
+        
+    def otf(self,x0,dx=0,dy=0,_caller='user'):
+        """
+        See __call__ for input arguments
+        Warning: result of otf will be unconsistent if undersampled!!!
+        This issue is solved with oversampling + binning in __call__ but not here
+        For the moment, the `_caller` keyword prevents user to misuse otf
+        """
+        
+        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)
+        
+        return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
+    
     def _otf_turbulent(self,x0):
         PSD,integral = self.psd(x0)
         L = self.system.D * self._samp_over
@@ -499,24 +510,18 @@ class ParametricPSFfromPSD(ParametricPSF):
     
     #@_lru_cache(maxsize=2)
     def _otf_diffraction(self):
-        Nx_over = self.npix[0]*self._k
-        Ny_over = self.npix[1]*self._k
-        
-        NpupX = _np.ceil(Nx_over/self._samp_over)
-        NpupY = _np.ceil(Ny_over/self._samp_over)
-        tab = _np.zeros((Nx_over, Ny_over), dtype=_np.complex)
+        nx, ny = self._nxnyOver
+        NpupX = _np.ceil(nx/self._samp_over)
+        NpupY = _np.ceil(ny/self._samp_over)
+        tab = _np.zeros((nx, ny), dtype=_np.complex)
         tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
         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):
-        Nx_over = self.npix[0]*self._k
-        Ny_over = self.npix[1]*self._k
-        
-        Y, X = _np.mgrid[0:Nx_over,0:Ny_over].astype(float)
-        X = (X-Nx_over/2) * 2*_np.pi*1j/Nx_over
-        Y = (Y-Ny_over/2) * 2*_np.pi*1j/Ny_over
-        return _np.exp(-X*dx*self._k - Y*dy*self._k)
+        Y, X = self._xyarray
+        ny,nx = X.shape
+        return _np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny))
     
     def __call__(self,x0,dx=0,dy=0):
         """