diff --git a/maoppy/example/create_muse_psf.py b/maoppy/example/create_muse_psf.py
new file mode 100644
index 0000000000000000000000000000000000000000..22abc527cdb4bc537c98b93f0b9bec8347ebd59f
--- /dev/null
+++ b/maoppy/example/create_muse_psf.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 24 15:18:38 2019
+
+@author: rfetick
+"""
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from maoppy.psfmodel import Psfao
+from maoppy.instrument import MUSE_NFM
+
+Npix = 128 # pixel size of PSF
+wvl = 600*1e-9 # wavelength [m]
+
+#%% 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)
+
+#%% Choose parameters and compute PSF
+r0 = 0.15 # Fried parameter [m]
+bck = 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,bck,amp,ax,beta]
+
+psf = Pmodel(param)
+
+#%% Plot results
+plt.figure(1)
+plt.clf()
+plt.pcolormesh(np.log(psf))
+plt.axis('image')
\ No newline at end of file