add jitter into ParametricPSFfromPSD

maoppy/psfmodel.py

@@ -202,6 +202,8 @@ class ParametricPSFfromPSD(ParametricPSF):
         self._fxy = None
         self._f2 = None
         self._otfDiffraction = None
+        self._jitter_mas = 0
+        self._otfJitter = 1
         
         self.fixed_k = fixed_k
         self.system = system
@@ -210,6 +212,12 @@ class ParametricPSFfromPSD(ParametricPSF):
         self._nparam = nparam
         self.bounds = ((-np.inf,)*nparam, (np.inf,)*nparam)
         
+      
+    def __str__(self):
+        s  = super().__str__ + '\n'
+        s += "jitter  : %.2f mas"%self.jitter_mas
+        return s  
+      
         
     @property
     def npix(self):
@@ -281,6 +289,24 @@ class ParametricPSFfromPSD(ParametricPSF):
         return fftshift(self._otfDiffraction)
     
     
+    @property
+    def jitter_mas(self):
+        return self._jitter_mas
+    
+    
+    @jitter_mas.setter
+    def jitter_mas(self,val):
+        self._jitter_mas = val
+        xx,yy = np.mgrid[0:self.npix[0],0:self.npix[1]]
+        dfx = 1/(self.npix[0]*self.system.resolution_mas)
+        dfy = 1/(self.npix[1]*self.system.resolution_mas)
+        xx = (xx - self.npix[0]//2)*dfx
+        yy = (yy - self.npix[1]//2)*dfy
+        # sigma = 1/val, but the writing below allows for val = 0
+        g = np.exp(-0.5*((xx*val)**2 + (yy*val)**2))
+        self._otfJitter = fftshift(g)
+    
+    
     def strehlOTF(self, parampsd):
         """Compute the Strehl ratio based on the sum of the OTF"""
         return np.real(np.sum(self._otf(parampsd))/np.sum(self._otfDiffraction))
@@ -368,7 +394,7 @@ class ParametricPSFfromPSD(ParametricPSF):
         else:
             otfTurb = self._otfTurbulent(parampsd, grad=False, **kwargs)
         
-        otf = otfTurb * self._otfDiffraction * otfShift
+        otf = otfTurb * self._otfDiffraction * otfShift * self._otfJitter
 
         if grad:
             return otf, gg