diff --git a/maoppy/data/muse_nfm.yml b/maoppy/data/muse_nfm.yml
index 821b03259fb90d89cd6a494ca3c4fdb67a4b39ba..c23e4aaf1099a0d78a720710578504544590f8cf 100644
--- a/maoppy/data/muse_nfm.yml
+++ b/maoppy/data/muse_nfm.yml
@@ -2,7 +2,7 @@ d: 8.0
 filters: {}
 fullname: VLT MUSE (NFM)
 gain: 5.0
-nact: 39
+nact: 34
 name: MUSE_NFM
 occ: 0.14
 phasemask_path: null
diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index 0a225005f439d26c9ddea1e6b05141acd4ea97df..e93435aa13695c4974ebf5909bd997470549e6f6 100644
--- a/maoppy/psfmodel.py
+++ b/maoppy/psfmodel.py
@@ -19,9 +19,21 @@ from maoppy.utils import binning as _binning
 _EPSILON = _np.sqrt(_sys.float_info.epsilon)
 
 #%% FUNCTIONS
-def oversample(samp):
-    """Find the minimal integer that allows oversampling"""
-    k = int(_np.ceil(2.0/samp))
+def oversample(samp, fixed_k = None):
+    """
+        Find the minimal integer that allows oversampling
+
+    Args:
+        samp (float): input sampling
+        fixed_k (int, optional): Oversampling factor to be fixed. Defaults to None.
+
+    Returns:
+        (k*samp, k): the oversampling (>=2) and the corresponding oversampling factor
+    """
+    if fixed_k == None: 
+        k = int(_np.ceil(2.0/samp))
+    else: 
+        k = fixed_k
     return (k*samp,k)
 
 #%% FITTING
@@ -396,17 +408,20 @@ class Gaussian(ParametricPSF):
 #%% MASTER CLASS
 class ParametricPSFfromPSD(ParametricPSF):
     """This class is NOT to be instantiated"""
-    def __init__(self,nparam,npix,system=None,samp=None):
+    def __init__(self,nparam,npix,system=None,samp=None, fixed_k = None):
         if not (type(npix) in [tuple,list,_np.ndarray]): raise ValueError("npix must be a tuple, list or numpy.ndarray")
         if len(npix)!=2: raise ValueError("npix must be of length = 2")
         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.fixed_k = fixed_k        
         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._nparam = nparam
         
+        
     @property
     def npix(self):
         return self._npix
@@ -421,9 +436,9 @@ class ParametricPSFfromPSD(ParametricPSF):
         return self._samp_over/self._k
     
     @samp.setter
-    def samp(self,value):
+    def samp(self, value):
         # Manage cases of undersampling
-        self._samp_over, self._k = oversample(value)
+        self._samp_over, self._k = oversample(value, fixed_k = self.fixed_k)
         self._computeXYarray()
         
     def _computeXYarray(self):
@@ -632,7 +647,7 @@ class Psfao(ParametricPSFfromPSD):
     ---------
     Fétick et al., August 2019, A&A, Vol.628
     """
-    def __init__(self,npix,system=None,Lext=10.,samp=None):
+    def __init__(self,npix,system=None,Lext=10.,samp=None, fixed_k = None):
         """
         Parameters
         ----------
@@ -646,7 +661,7 @@ class Psfao(ParametricPSFfromPSD):
             Von-Karman external scale (default = 10 m)
             Useless if Fao >> 1/Lext
         """
-        super().__init__(7,npix,system=system,samp=samp)
+        super().__init__(7,npix,system=system,samp=samp,fixed_k=fixed_k)
         
         self.Lext = Lext
         # r0,C,A,alpha,ratio,theta,beta 
@@ -688,6 +703,7 @@ class Psfao(ParametricPSFfromPSD):
         
         moff = moffat(f2D,param,norm=Fao) * (F2 < Fao**2.)
         moff[nx0,ny0] = 0.0 # Set Moffat PSD = 0 at null frequency
+        # Newly added for the PSFAO19 model
         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