From 4257749292c229531a0673e17b2062bbda17c915 Mon Sep 17 00:00:00 2001
From: Romain Fetick <rfetick@OSU.PYTHEAS>
Date: Thu, 11 Jul 2024 13:57:17 +0200
Subject: [PATCH] update PsfaoAliasing

---
 maoppy/psfmodel.py | 41 +++++++++++++++++++++++++++++++++--------
 1 file changed, 33 insertions(+), 8 deletions(-)

diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index ca4345d..7c19ac6 100644
--- a/maoppy/psfmodel.py
+++ b/maoppy/psfmodel.py
@@ -653,7 +653,7 @@ class Psfao(ParametricPSFfromPSD):
             
         Keywords
         --------
-        system : OpticalSystem
+        system : Instrument
             Optical system for this PSF
         samp : float
             Sampling at the observation wavelength
@@ -851,18 +851,43 @@ class Psfao(ParametricPSFfromPSD):
 
 
 class PsfaoAliasing(Psfao):
-    """Ignore background and use aliasing instead"""
+    """
+    Same as Psfao but without background parameter.
+    The background is automatically set using `r0` and `system.aliasing`.
+    
+    Note
+    ----
+    The `bounds` property holds 7 parameters since it is called by the parent
+    class `Psfao.psd(...)` with 7 parameters and not by `PsfaoAliasing.psd(...)`.
+    """
+        
+    @property
+    def param_name(self):
+        pn = super().param_name
+        return [pn[0]] + pn[2:]
+    
+    @property
+    def param_comment(self):
+        pc = super().param_comment
+        return [pc[0]] + pc[2:]
     
     def psd(self, parampsd, *args, grad=False, **kwargs):
-        # parampsd can be called with a dictionary
+        """
+        Same function as `Psfao.psd(...)`.
+        `parampsd` only holds 6 parameters, the background one being removed.
+        """
+        # parampsd can be called with a dictionary (6 parameters)
         if isinstance(parampsd,dict):
             parampsd = self.dict2list(parampsd)
         
         r0 = parampsd[0]
-        bck = self.system.aliasing*0.0229*6/5* r0**(-5/3) * self.cutoffAOfreq**(-11/3)
-        parampsd[1] = bck
+        var_fitting = 0.023*6*np.pi/5 * (r0*self.cutoffAOfreq)**(-5./3.)
+        var_alias = self.system.aliasing * var_fitting # [rad²]
+        surf_corr = np.pi * self.cutoffAOfreq**2 # [1/m²]
+        bck = var_alias / surf_corr
+        parampsd_seven = np.concatenate(([r0,bck],parampsd[1:]))
         
-        out = super().psd(parampsd, *args, grad=grad, **kwargs)
+        out = super().psd(parampsd_seven, *args, grad=grad, **kwargs)
         
         if grad:
             # I need to update the gradient and its integral!
@@ -871,8 +896,8 @@ class PsfaoAliasing(Psfao):
             newbckmsk = dbck_dr0*self._maskin
             gg[0,...] += newbckmsk
             integral_g[0] += np.sum(newbckmsk) * self.pix2freq**2
-            gg[1,...] = 0
-            integral_g[1] = 0
+            gg = np.delete(gg, (1), axis=0)
+            integral_g = np.delete(integral_g, (1), axis=0)
             return psd, integral, gg, integral_g
                 
         return out
-- 
GitLab