Commit f56044d5 authored by FETICK Romain's avatar FETICK Romain
Browse files

Update and add one example

parent 97640621
......@@ -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
"""
......
#!/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')
......@@ -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
"""
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment