add TurbulenceJitter class

maoppy/psfmodel.py

@@ -724,6 +724,96 @@ class Turbulent(ParametricPSFfromPSD):
         return PSD, integral
 
 
+class TurbulentJitter(ParametricPSFfromPSD):
+    """
+    Summary
+    -------
+    Same as `Turbulent` but with extra gaussian jitter on one axis
+    
+    Warning
+    -------
+    Not finished yet, do not use!
+    #TODO: finish to code the class
+    """
+    def __init__(self,npix,system=None,samp=None):
+        """
+        Parameters
+        ----------
+        npix : tuple
+            Size of output PSF
+        system : OpticalSystem
+            Optical system for this PSF
+        samp : float
+            Sampling at the observation wavelength
+        """
+        print("Warning: TurbulentJitter is not finished yet, do not use!")
+        super().__init__(4,npix,system=system,samp=samp)
+        
+        # r0,L0,sigma,theta
+        bounds_down = [_EPSILON, _EPSILON, _EPSILON, -_np.inf]
+        bounds_up = [_np.inf, _np.inf, _np.inf, _np.inf]
+        self.bounds = (bounds_down,bounds_up)
+        
+    def psd(self,x0,grad=False):
+        """Compute the PSD model from parameters
+        PSD is given in [rad²/f²] = [rad² m²]
+        
+        Parameters
+        ----------
+        x0 : numpy.array (dim=1), tuple, list
+            See __doc__ for more details
+            
+        Returns
+        -------
+        psd : numpy.array (dim=2)
+        integral : float : the integral of the `psd` up to infinity
+            
+        """
+        self.check_parameters(x0)
+        nx0,ny0 = self._nullFreqIndex
+        f2D = self._fxfy
+        
+        F2 = f2D[0]**2. + f2D[1]**2.
+        
+        r0,Lext = x0
+        PSD = 0.0229* r0**(-5./3.) * ((1./Lext**2.) + F2)**(-11./6.)
+
+        # Set PSD = 0 at null frequency
+        PSD[nx0,ny0] = 0.0
+        
+        # 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 (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
+        
+        if grad:
+            g = _np.zeros((len(x0),)+F2.shape)
+            g[0,...] = PSD*(-5./3)/r0
+            g[1,...] = PSD*(-11./6)/((1./Lext**2.) + F2)*(-2/Lext**3)
+            
+            # compute integral gradient
+            integral_g = [0,0]
+            for i in range(2): integral_g[i] = _np.sum(g[i,...]*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum (in the array's tangent circle)
+            integral_g[0] += integral_out*(-5./3)/r0
+        
+        if grad: return PSD,integral,g,integral_g
+        return PSD, integral
+    
+    def otf(self,x0,dx=0,dy=0,_caller='user',grad=False):
+        if grad:
+            raise NotImplementedError("gradient computation not implemented yet")
+        else:
+            otfsuper = super().otf(x0[0:2],dx=dx,dy=dy,_caller=_caller,grad=grad)
+            return otfsuper * self._otfJitter(x0[2:])
+    
+    def _otfJitter(self,x0):
+        """x0 = [sigma,theta]"""
+        ny,nx = self._xyarray[0].shape
+        x = _np.arange(nx)-nx/2
+        u = _np.zeros((ny,nx))
+        
+
 #%% PSFAO MODEL        
 class Psfao(ParametricPSFfromPSD):
     """