Add Gaussian PSF

maoppy/example/fit_muse_simulated_psf.py

@@ -19,7 +19,7 @@ from maoppy.instrument import MUSE_NFM
 npix = 128 # pixel size of PSF
 wvl = 600*1e-9 # wavelength [m]
 
-flux = 1e6
+flux = 1e5
 ron = MUSE_NFM.ron
 bckgd = 1.0
 
@@ -38,7 +38,8 @@ param = [r0,b,amp,ax,beta]
 
 psf = Pmodel(param,dx=0,dy=0)
 
-image = flux*psf + bckgd + ron*np.random.randn(npix,npix)
+noise = ron*np.random.randn(npix,npix)
+image = flux*psf + bckgd + noise
 
 #%% Fit the PSF
 
@@ -59,7 +60,8 @@ print("%9s %10.4g %10.4g"%('Bck [e-]',bckgd,bck_fit))
 print("%9s %10.2f %10.2f"%('r0 [cm]',r0*100,out.x[0]*100))
 print("%9s %10.2f %10.2f"%('A [rad²]',amp,out.x[2]))
 print(31*'-')
-print("shape RMS error = %.2f%%"%(100*rmserror(fitao,image)))
+print("shape + noise RMS error = %.2f%%"%(100*rmserror(fitao,image)))
+print("        noise RMS error = %.2f%%"%(np.sqrt(np.sum(noise**2))/np.sum(image)*100))
 print(31*'-')
 
 vmin,vmax = np.arcsinh([image.min(),image.max()])

maoppy/psfmodel.py

@@ -76,6 +76,7 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
         Initial guess for parameters
     weights : numpy.ndarray
         Least-square weighting matrix (same size as `psf`)
+        Inverse of the noise's variance
         Default: uniform weighting
     dxdy : tuple of two floats
         Eventual guess on PSF shifting
@@ -354,6 +355,31 @@ class Moffat(ParametricPSF):
         super(Moffat,self).tofits(param,filename,*args,keys=keys,**kwargs)
 
 
+class Gaussian(ParametricPSF):
+    def __init__(self,Npix):
+        self.Npix = Npix
+        bounds_down = [_EPSILON,_EPSILON,-np.inf]
+        bounds_up = [np.inf for i in range(3)]
+        self.bounds = (bounds_down,bounds_up)
+        
+        self._XY = np.mgrid[0:Npix[0],0:Npix[1]]
+        self._XY[1] -= Npix[0]/2
+        self._XY[0] -= Npix[1]/2
+    
+    def __call__(self,x,dx=0,dy=0):
+        """
+        Parameters
+        ----------
+        x : list, tuple, numpy.ndarray (len=4)
+            x[0] - Sigma X
+            x[1] - Sigma Y
+            x[2] - Theta
+        """
+        y = np.concatenate((x,[dx,dy]))
+        return gauss(self._XY,y)
+    
+    def tofits(self,param,filename,*args,keys=["SIGMA_X","SIGMA_Y","THETA"],**kwargs):
+        super(Moffat,self).tofits(param,filename,*args,keys=keys,**kwargs)
 
 
 def oversample(samp):