compute Psfao derivatives, remove numerical normalization

maoppy/example/psfao_derivative.py

@@ -25,7 +25,7 @@ Pmodel = Psfao((Npix,Npix),system=muse_nfm,samp=samp)
 r0 = 0.15 # Fried parameter [m]
 bck = 1e-5 # Phase PSD background [rad² m²]
 amp = 5.0 # Phase PSD Moffat amplitude [rad²]
-alpha = 1 # Phase PSD Moffat alpha [1/m]
+alpha = 0.8 # Phase PSD Moffat alpha [1/m]
 ratio = 1.2
 theta = np.pi/4
 beta = 1.6 # Phase PSD Moffat beta power law

maoppy/psfmodel.py

@@ -282,7 +282,7 @@ def moffat(XY, param, norm=None, removeInside=0, grad=False):
         Fout = (1 +         (norm**2)/(ax*ay))**(1-beta)
         Fin  = (1 + (removeInside**2)/(ax*ay))**(1-beta)
         F = 1/(Fin-Fout)
-        if grad: # TODO: debug gradient for norm=float
+        if grad:
             dE = [-E/ax,-E/ay,0,E/(beta-1)]
             k = (1-beta)*Fout/(1 + (norm**2)/(ax*ay))*(norm**2)/(ax*ay)
             dFout = _np.array([-k/ax,-k/ay,0,-Fout*_np.log(1 + (norm**2)/(ax*ay))]) if norm<_np.inf else _np.zeros(4)
@@ -772,6 +772,7 @@ class Psfao(ParametricPSFfromPSD):
         self.check_parameters(x0)
         nx0,ny0 = self._nullFreqIndex
         f2D = self._fxfy
+        pix = self._pix2freq
         F2 = f2D[0]**2. + f2D[1]**2.
         Fao = self.system.Nact/(2.0*self.system.D)
         maskin = (F2 < Fao**2.)
@@ -785,18 +786,24 @@ class Psfao(ParametricPSFfromPSD):
         ay = alpha/ratio
         param = (ax,ay,theta,beta,0,0)
         
-        removeInside = (1+_np.sqrt(2))/2 * self._pix2freq/2 # remove central pixel in energy computation
-        moff = moffat(f2D,param,norm=Fao,removeInside=removeInside) * maskin
+        removeInside = (1+_np.sqrt(2))/2 * pix/2 # remove central pixel in energy computation
+        if grad:
+            moff,dm = moffat(f2D,param,norm=Fao,removeInside=removeInside,grad=True)
+        else:
+            moff = moffat(f2D,param,norm=Fao,removeInside=removeInside)
+        moff *= maskin
+        
+        # Warning: Numerical normalization! (not compatible with analytical gradient)
         moff[nx0,ny0] = 0.0 # Set Moffat PSD = 0 at null frequency
-        moff = moff / (_np.sum(moff)*self._pix2freq**2)  # normalize moffat numerically to get correct A=sigma² in the AO zone
+        #moff = moff / (_np.sum(moff)*pix**2)  # normalize moffat numerically to get correct A=sigma² in the AO zone
         # Warning: Moffat numerical normalization generates strehlOTF jump when "_k" is changed
         
         PSD += maskin * (C + A*moff) # before I wrote "moffat(f2D,param,norm=Fao)" instead of "moff", so no numerical normalization
         PSD[nx0,ny0] = 0.0 # Set PSD = 0 at null frequency
         
         # 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
+        fmax = _np.min([nx0,ny0])*pix
+        integral_in = _np.sum(PSD*(F2<(fmax**2))) * pix**2 # numerical sum
         integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
         integral = integral_in + integral_out
         
@@ -810,9 +817,13 @@ class Psfao(ParametricPSFfromPSD):
             # derivative towards A
             g[2,...] = maskin*moff
             # derivative towards alpha
+            g[3,...] = A*maskin*(dm[0,...]*ratio + dm[1,...]/ratio)
             # derivative towards ratio
+            g[4,...] = A*maskin*(dm[0,...]*alpha - dm[1,...]*ay/ratio)
             # derivative towards theta
+            g[5,...] = A*maskin*dm[2,...]
             # derivative towards beta
+            g[6,...] = A*maskin*dm[3,...]
             # Remove central freq from all derivative
             g[:,nx0,ny0] = 0