From f56044d53b7c7615be7d2eaf0bf935ed90865210 Mon Sep 17 00:00:00 2001
From: rfetick <r.fetick@gmail.com>
Date: Fri, 8 May 2020 19:18:30 +0200
Subject: [PATCH] Update and add one example

---
 maoppy/example/create_muse_psf.py        |  2 +
 maoppy/example/fit_muse_simulated_psf.py | 84 ++++++++++++++++++++++++
 maoppy/example/fit_zimpol_psf.py         |  3 +
 maoppy/psfmodel.py                       |  6 +-
 4 files changed, 92 insertions(+), 3 deletions(-)
 create mode 100644 maoppy/example/fit_muse_simulated_psf.py

diff --git a/maoppy/example/create_muse_psf.py b/maoppy/example/create_muse_psf.py
index f74e1d1..2e886f1 100644
--- a/maoppy/example/create_muse_psf.py
+++ b/maoppy/example/create_muse_psf.py
@@ -3,6 +3,8 @@
 """
 Created on Wed Jul 24 15:18:38 2019
 
+Generate a VLT/MUSE-NFM PSF with the Psfao model.
+
 @author: rfetick
 """
 
diff --git a/maoppy/example/fit_muse_simulated_psf.py b/maoppy/example/fit_muse_simulated_psf.py
new file mode 100644
index 0000000..1e92f2d
--- /dev/null
+++ b/maoppy/example/fit_muse_simulated_psf.py
@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May  8 18:45:28 2020
+
+Generate an observation of a star with MUSE-NFM, at a given wavelength
+Fit the simulated image with the Psfao model
+Plot results and display fitting accuracy
+
+@author: rfetick
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from maoppy.psfmodel import Psfao, psffit, rmserror
+from maoppy.instrument import MUSE_NFM
+
+npix = 128 # pixel size of PSF
+wvl = 600*1e-9 # wavelength [m]
+
+flux = 1e6
+ron = MUSE_NFM.ron
+bckgd = 1.0
+
+#%% Initialize PSF model
+samp = MUSE_NFM.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
+Pmodel = Psfao((npix,npix),system=MUSE_NFM,symmetric=True,samp=samp)
+
+#%% Generate a MUSE-NFM PSF
+r0 = 0.15 # Fried parameter [m]
+b = 1e-7 # Phase PSD background [rad² m²]
+amp = 1.4 # Phase PSD Moffat amplitude [rad²]
+ax = 0.1 # Phase PSD Moffat alpha [1/m]
+beta = 1.6 # Phase PSD Moffat beta power law
+
+param = [r0,b,amp,ax,beta]
+
+psf = Pmodel(param,dx=0,dy=0)
+
+image = flux*psf + bckgd + ron*np.random.randn(npix,npix)
+
+#%% Fit the PSF
+
+guess = [0.145,2e-7,1.2,0.08,1.5]
+w = np.ones_like(image)/ron**2.0
+
+out = psffit(image,Psfao,guess,weights=w,system=MUSE_NFM,samp=samp,symmetric=True)
+
+flux_fit, bck_fit = out.flux_bck
+fitao = flux_fit*out.psf + bck_fit
+
+#%% Plots
+print(31*'-')
+print("%9s %10s %10s"%('','MUSE-NFM','Psfao'))
+print(31*'-')
+print("%9s %10.4g %10.4g"%('Flux [e-]',flux,flux_fit))
+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(31*'-')
+
+vmin,vmax = np.arcsinh([image.min(),image.max()])
+
+plt.figure(1)
+plt.clf()
+
+plt.subplot(131)
+plt.imshow(np.arcsinh(image),vmin=vmin,vmax=vmax)
+plt.axis('off')
+plt.title('MUSE-NFM simu')
+
+plt.subplot(132)
+plt.imshow(np.arcsinh(fitao),vmin=vmin,vmax=vmax)
+plt.axis('off')
+plt.title('Psfao fit')
+
+plt.subplot(133)
+plt.imshow(np.arcsinh(image-fitao))
+plt.axis('off')
+plt.title('residuals')
+
diff --git a/maoppy/example/fit_zimpol_psf.py b/maoppy/example/fit_zimpol_psf.py
index ec74172..706190f 100644
--- a/maoppy/example/fit_zimpol_psf.py
+++ b/maoppy/example/fit_zimpol_psf.py
@@ -3,6 +3,9 @@
 """
 Created on Wed May 29 11:06:34 2019
 
+Fit a VLT/SPHERE/Zimpol PSF with the Psfao model.
+You need a Zimpol PSF in FITS format beforehand.
+
 @author: rfetick
 """
 
diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index 4cd98df..ee48c2a 100644
--- a/maoppy/psfmodel.py
+++ b/maoppy/psfmodel.py
@@ -12,15 +12,15 @@ from astropy.io import fits
 import time
 from numpy.fft import fft2, fftshift, ifft2
 from maoppy.config import _EPSILON
-from maoppy.utils import binning, getgaussnoise
+from maoppy.utils import binning #, getgaussnoise
 
 #%% FITTING FUNCTION
 def rmserror(model,image,weights=1.0,background=True, positive_bck=False):
     """Compute the RMS error between a PSF model and an image"""
     amp,bck = lsq_flux_bck(model, image, weights=weights, background=background, positive_bck=positive_bck)
     diff = amp*model+bck - image
-    noise_std = getgaussnoise(image)
-    err = np.sqrt(np.sum(diff**2.0 - noise_std**2.0))
+    #noise_std = getgaussnoise(image)
+    err = np.sqrt(np.sum(diff**2.0))
     err = err/np.sum(image)
     return err
 
-- 
GitLab