From 2ac02eb3acecacd1155ce8b0835cf51442ba156b Mon Sep 17 00:00:00 2001
From: rfetick <r.fetick@gmail.com>
Date: Fri, 8 May 2020 20:11:26 +0200
Subject: [PATCH] Add Gaussian PSF
---
maoppy/example/fit_muse_simulated_psf.py | 8 +++++---
maoppy/psfmodel.py | 26 ++++++++++++++++++++++++
2 files changed, 31 insertions(+), 3 deletions(-)
diff --git a/maoppy/example/fit_muse_simulated_psf.py b/maoppy/example/fit_muse_simulated_psf.py
index 1e92f2d..6686d66 100644
--- a/maoppy/example/fit_muse_simulated_psf.py
+++ b/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()])
diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index ee48c2a..7ac1f02 100644
--- a/maoppy/psfmodel.py
+++ b/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):
--
GitLab