Skip to content
Snippets Groups Projects
Commit e5bc3070 authored by rfetick's avatar rfetick
Browse files

Dynamically defined Instruments as attributes

parent 3c9cc925
No related branches found
No related tags found
No related merge requests found
...@@ -13,16 +13,14 @@ import matplotlib.pyplot as plt ...@@ -13,16 +13,14 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
from maoppy.psfmodel import Psfao from maoppy.psfmodel import Psfao
from maoppy.instrument import load from maoppy.instrument import muse_nfm
Npix = 128 # pixel size of PSF Npix = 128 # pixel size of PSF
wvl = 600*1e-9 # wavelength [m] wvl = 600*1e-9 # wavelength [m]
MUSE_NFM = load('muse_nfm')
#%% Initialize PSF model #%% Initialize PSF model
samp = MUSE_NFM.samp(wvl) # sampling (2.0 for Shannon-Nyquist) samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((Npix,Npix),system=MUSE_NFM,samp=samp) Pmodel = Psfao((Npix,Npix),system=muse_nfm,samp=samp)
#%% Choose parameters and compute PSF #%% Choose parameters and compute PSF
r0 = 0.15 # Fried parameter [m] r0 = 0.15 # Fried parameter [m]
......
...@@ -18,20 +18,18 @@ import matplotlib.pyplot as plt ...@@ -18,20 +18,18 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
from maoppy.psfmodel import Psfao, psffit, rmserror from maoppy.psfmodel import Psfao, psffit, rmserror
from maoppy.instrument import load from maoppy.instrument import muse_nfm
MUSE_NFM = load('muse_nfm')
npix = 128 # pixel size of PSF npix = 128 # pixel size of PSF
wvl = 600*1e-9 # wavelength [m] wvl = 600*1e-9 # wavelength [m]
flux = 1e5 flux = 1e5
ron = MUSE_NFM.ron ron = muse_nfm.ron
bckgd = 1.0 bckgd = 1.0
#%% Initialize PSF model #%% Initialize PSF model
samp = MUSE_NFM.samp(wvl) # sampling (2.0 for Shannon-Nyquist) samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((npix,npix),system=MUSE_NFM,samp=samp) Pmodel = Psfao((npix,npix),system=muse_nfm,samp=samp)
#%% Generate a MUSE-NFM PSF #%% Generate a MUSE-NFM PSF
r0 = 0.15 # Fried parameter [m] r0 = 0.15 # Fried parameter [m]
...@@ -54,7 +52,7 @@ guess = [0.145,2e-7,1.2,0.08,ratio,theta,1.5] ...@@ -54,7 +52,7 @@ guess = [0.145,2e-7,1.2,0.08,ratio,theta,1.5]
w = np.ones_like(image)/ron**2.0 w = np.ones_like(image)/ron**2.0
fixed = [False,False,False,False,True,True,False] fixed = [False,False,False,False,True,True,False]
out = psffit(image,Psfao,guess,weights=w,system=MUSE_NFM,samp=samp,fixed=fixed) out = psffit(image,Psfao,guess,weights=w,system=muse_nfm,samp=samp,fixed=fixed)
flux_fit, bck_fit = out.flux_bck flux_fit, bck_fit = out.flux_bck
fitao = flux_fit*out.psf + bck_fit fitao = flux_fit*out.psf + bck_fit
......
...@@ -8,7 +8,7 @@ Show the possibility of fitting a PSF with Psfao on a small field of view ...@@ -8,7 +8,7 @@ Show the possibility of fitting a PSF with Psfao on a small field of view
and to retrieve the PSF shape on an increased field of view. and to retrieve the PSF shape on an increased field of view.
-------------- --------------
This script is a more complex version of 'fit_muse_simulated_psf.py' This script is a more complex version of 'fit_muse_simulated_psf.py'
The other one should be run first. The other one should be run first for better understanding.
-------------- --------------
Generate an observation of a star with MUSE-NFM, at a given wavelength Generate an observation of a star with MUSE-NFM, at a given wavelength
...@@ -26,20 +26,18 @@ from matplotlib.patches import Rectangle ...@@ -26,20 +26,18 @@ from matplotlib.patches import Rectangle
import numpy as np import numpy as np
from maoppy.psfmodel import Psfao, psffit, rmserror from maoppy.psfmodel import Psfao, psffit, rmserror
from maoppy.instrument import load from maoppy.instrument import muse_nfm
MUSE_NFM = load('muse_nfm')
npix = 128 # pixel size of PSF npix = 128 # pixel size of PSF
wvl = 600*1e-9 # wavelength [m] wvl = 600*1e-9 # wavelength [m]
flux = 1e6 flux = 1e6
ron = MUSE_NFM.ron ron = muse_nfm.ron
bckgd = 1.0 bckgd = 1.0
#%% Initialize PSF model #%% Initialize PSF model
samp = MUSE_NFM.samp(wvl) # sampling (2.0 for Shannon-Nyquist) samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((npix,npix),system=MUSE_NFM,samp=samp) Pmodel = Psfao((npix,npix),system=muse_nfm,samp=samp)
#%% Generate a MUSE-NFM PSF #%% Generate a MUSE-NFM PSF
r0 = 0.15 # Fried parameter [m] r0 = 0.15 # Fried parameter [m]
...@@ -68,7 +66,7 @@ w = np.ones_like(im_crop)/ron**2.0 ...@@ -68,7 +66,7 @@ w = np.ones_like(im_crop)/ron**2.0
fixed = [False,False,False,False,True,True,False] fixed = [False,False,False,False,True,True,False]
out = psffit(im_crop,Psfao,guess,weights=w,fixed=fixed,npixfit=npix, # fit keywords out = psffit(im_crop,Psfao,guess,weights=w,fixed=fixed,npixfit=npix, # fit keywords
system=MUSE_NFM,samp=samp) # Psfao keywords system=muse_nfm,samp=samp) # Psfao keywords
flux_fit, bck_fit = out.flux_bck flux_fit, bck_fit = out.flux_bck
fitao = flux_fit*out.psf + bck_fit fitao = flux_fit*out.psf + bck_fit
......
...@@ -15,13 +15,11 @@ import numpy as np ...@@ -15,13 +15,11 @@ import numpy as np
from maoppy.utils import imcenter, circavgplt from maoppy.utils import imcenter, circavgplt
from maoppy.psfmodel import psffit, Psfao from maoppy.psfmodel import psffit, Psfao
from maoppy.instrument import load from maoppy.instrument import zimpol
ZIMPOL = load('zimpol')
# %% PARAMETERS TO MODIFY # %% PARAMETERS TO MODIFY
psf_path = "path/to/your/zimpol_psf.fits" # set the path to your PSF *.fits file psf_path = "path/to/your/zimpol_psf.fits" # set the path to your PSF *.fits file
filt = ZIMPOL.filters["N_R"] # ZIMPOL filter filt = zimpol.filters["N_R"] # ZIMPOL filter
Npix = 512 # size of PSF for fitting [pixels]. The larger the better Npix = 512 # size of PSF for fitting [pixels]. The larger the better
r0 = 0.18 # Fried parameter [m] r0 = 0.18 # Fried parameter [m]
...@@ -37,11 +35,11 @@ fitsPSF = fits.open(psf_path) ...@@ -37,11 +35,11 @@ fitsPSF = fits.open(psf_path)
hdr = fitsPSF[0].header hdr = fitsPSF[0].header
psf = fitsPSF[0].data psf = fitsPSF[0].data
psf = imcenter(psf, (Npix, Npix), maxi=True) * ZIMPOL.gain psf = imcenter(psf, (Npix, Npix), maxi=True) * zimpol.gain
# %% FIT PSFAO model # %% FIT PSFAO model
wvl = filt[0] wvl = filt[0]
samp = ZIMPOL.samp(wvl) samp = zimpol.samp(wvl)
x0 = [r0,moff_C,moff_A,moff_alpha,moff_ratio,moff_theta,moff_beta] x0 = [r0,moff_C,moff_A,moff_alpha,moff_ratio,moff_theta,moff_beta]
fixed = [False,False,False,False,True,True,False] fixed = [False,False,False,False,True,True,False]
...@@ -49,7 +47,7 @@ fixed = [False,False,False,False,True,True,False] ...@@ -49,7 +47,7 @@ fixed = [False,False,False,False,True,True,False]
weights = np.ones_like(psf) weights = np.ones_like(psf)
print("PSFAO fitting (please wait)") print("PSFAO fitting (please wait)")
out = psffit(psf,Psfao,x0,weights=weights,system=ZIMPOL,samp=samp,fixed=fixed) out = psffit(psf,Psfao,x0,weights=weights,system=zimpol,samp=samp,fixed=fixed)
out.x[5] = out.x[5]%np.pi # make angle between 0 and PI out.x[5] = out.x[5]%np.pi # make angle between 0 and PI
......
...@@ -6,31 +6,34 @@ Created on Mon May 27 17:30:51 2019 ...@@ -6,31 +6,34 @@ Created on Mon May 27 17:30:51 2019
@author: rfetick @author: rfetick
""" """
from maoppy.utils import circarr, RAD2ARCSEC from maoppy.utils import circarr as _circarr
import yaml from maoppy.utils import RAD2ARCSEC as _RAD2ARCSEC
import os
from astropy.io import fits import sys as _sys
import numpy as np import yaml as _yaml
from scipy.interpolate import interp2d import os as _os
from astropy.io import fits as _fits
import numpy as _np
from scipy.interpolate import interp2d as _interp2d
def _get_data_folder(): def _get_data_folder():
folder = os.path.abspath(__file__) folder = _os.path.abspath(__file__)
folder = os.sep.join(folder.split(os.sep)[0:-1])+os.sep+'data'+os.sep folder = _os.sep.join(folder.split(_os.sep)[0:-1])+_os.sep+'data'+_os.sep
return folder return folder
def load(name): #def load(name):
"""Load an instrument saved in the `data` folder of MAOPPY""" # """Load an instrument saved in the `data` folder of MAOPPY"""
folder = _get_data_folder() # folder = _get_data_folder()
with open(folder+name.lower()+".yml",'r') as f: # with open(folder+name.lower()+".yml",'r') as f:
d = yaml.full_load(f) # d = _yaml.full_load(f)
instru = Instrument(D=d['d'],occ=d['occ'],res=d['res'],gain=d['gain'],ron=d['ron'],Nact=d['nact']) # instru = Instrument(D=d['d'],occ=d['occ'],res=d['res'],gain=d['gain'],ron=d['ron'],Nact=d['nact'])
instru.name = d['name'] # instru.name = d['name']
instru.fullname = d['fullname'] # instru.fullname = d['fullname']
instru.filters = d['filters'] # instru.filters = d['filters']
instru.phasemask_path = d['phasemask_path'] # instru.phasemask_path = d['phasemask_path']
instru.phasemask_shift = d['phasemask_shift'] # instru.phasemask_shift = d['phasemask_shift']
return instru # return instru
#%% INSTRUMENT CLASS #%% INSTRUMENT CLASS
class Instrument(object): class Instrument(object):
...@@ -106,25 +109,25 @@ class Instrument(object): ...@@ -106,25 +109,25 @@ class Instrument(object):
@property @property
def resolution_mas(self): def resolution_mas(self):
return self.resolution_rad * RAD2ARCSEC * 1e3 return self.resolution_rad * _RAD2ARCSEC * 1e3
def pupil(self,Npix,wvl=None,samp=None): def pupil(self,Npix,wvl=None,samp=None):
"""Returns the 2D array of the pupil transmission function (complex data)""" """Returns the 2D array of the pupil transmission function (complex data)"""
Dpix = min(Npix)/2 Dpix = min(Npix)/2
pup = circarr(Npix) pup = _circarr(Npix)
if self.phasemask_enable: if self.phasemask_enable:
if self._phasemask is None: if self._phasemask is None:
if self.phasemask_path is None: raise ValueError('phasemask_path must be defined') if self.phasemask_path is None: raise ValueError('phasemask_path must be defined')
p = fits.open(self.phasemask_path)[0].data * 1e-9 # fits data in nm, converted here to meter p = _fits.open(self.phasemask_path)[0].data * 1e-9 # fits data in nm, converted here to meter
x = np.arange(p.shape[0])/p.shape[0] x = _np.arange(p.shape[0])/p.shape[0]
y = np.arange(p.shape[1])/p.shape[1] y = _np.arange(p.shape[1])/p.shape[1]
self._phasemask = interp2d(x,y,p) self._phasemask = _interp2d(x,y,p)
cx,cy = self.phasemask_shift cx,cy = self.phasemask_shift
x = np.arange(Npix[0])/Npix[0] - cx/Npix[0] x = _np.arange(Npix[0])/Npix[0] - cx/Npix[0]
y = np.arange(Npix[1])/Npix[1] - cy/Npix[1] y = _np.arange(Npix[1])/Npix[1] - cy/Npix[1]
wf = self._phasemask(x,y) wf = self._phasemask(x,y)
if wvl is None: wvl = self.wvl(samp) # samp must be defined if wvl is None if wvl is None: wvl = self.wvl(samp) # samp must be defined if wvl is None
wf = np.exp(2j*np.pi/wvl*wf) wf = _np.exp(2j*_np.pi/wvl*wf)
else: else:
wf = 1.0 + 0j # complex type for output array, even if real data wf = 1.0 + 0j # complex type for output array, even if real data
return (pup < Dpix) * (pup >= Dpix*self.occ) * wf return (pup < Dpix) * (pup >= Dpix*self.occ) * wf
...@@ -151,18 +154,23 @@ class Instrument(object): ...@@ -151,18 +154,23 @@ class Instrument(object):
data['filters'] = self.filters data['filters'] = self.filters
data['phasemask_path'] = self.phasemask_path data['phasemask_path'] = self.phasemask_path
data['phasemask_shift'] = self.phasemask_shift data['phasemask_shift'] = self.phasemask_shift
with open(filename,'w') as f: yaml.dump(data,f) with open(filename,'w') as f: _yaml.dump(data,f)
#%% OLD WAY to LOAD INSTRUMENTS (now use the `load` function) #%% LOAD INSTRUMENT INSTANCES (make them attributes of this module)
_this_module = _sys.modules[__name__]
_d = _get_data_folder()
_l = [_f for _f in _os.listdir(_d) if _f.endswith('.yml')]
for _f in _l:
with open(_d+_f,'r') as _fi:
_i = _yaml.full_load(_fi)
_instru = Instrument(D=_i['d'],occ=_i['occ'],res=_i['res'],gain=_i['gain'],ron=_i['ron'],Nact=_i['nact'])
_instru.name = _i['name']
_instru.fullname = _i['fullname']
_instru.filters = _i['filters']
_instru.phasemask_path = _i['phasemask_path']
_instru.phasemask_shift = _i['phasemask_shift']
_n = _i['name'].lower().replace(" ","_") # format name
setattr(_this_module,_n,_instru)
del _this_module, _d, _l, _f, _instru, _fi, _i, _n
#ZIMPOL = Instrument(D=8.,occ=0.14,res=30*1e-6/1768.,gain=10.5,ron=20.,Nact=40)
#ZIMPOL.filters["V"] = (554*1e-9, 80.6*1e-9)
#ZIMPOL.filters["N_R"] = (645.9*1e-9, 56.7*1e-9)
#ZIMPOL.name = "ZIMPOL"
#ZIMPOL.fullname = "VLT SPHERE/ZIMPOL"
#
#
#MUSE_NFM = Instrument(D=8.,occ=0.14,res=237.15*1e-6/1980.,gain=5.,ron=15.,Nact=39)
#MUSE_NFM.name = "MUSE_NFM"
#MUSE_NFM.fullname = "VLT MUSE (NFM)"
...@@ -9,7 +9,7 @@ Created on Fri May 8 20:21:20 2020 ...@@ -9,7 +9,7 @@ Created on Fri May 8 20:21:20 2020
import unittest import unittest
import numpy as np import numpy as np
from maoppy.psfmodel import lsq_flux_bck, moffat, gauss, Psfao from maoppy.psfmodel import lsq_flux_bck, moffat, gauss, Psfao
from maoppy.instrument import load from maoppy.instrument import zimpol, muse_nfm
class TestLsqFluxBck(unittest.TestCase): class TestLsqFluxBck(unittest.TestCase):
def test_fluxbck(self): def test_fluxbck(self):
...@@ -55,7 +55,7 @@ class TestPsfao(unittest.TestCase): ...@@ -55,7 +55,7 @@ class TestPsfao(unittest.TestCase):
def test_oversampling(self): def test_oversampling(self):
npix = 100 npix = 100
samp = 0.3 samp = 0.3
P = Psfao((npix,npix),system=load('muse_nfm'),samp=samp) P = Psfao((npix,npix),system=muse_nfm,samp=samp)
self.assertEqual(samp,P.samp) self.assertEqual(samp,P.samp)
self.assertGreaterEqual(P._samp_over,2.0) self.assertGreaterEqual(P._samp_over,2.0)
self.assertEqual(P._k%1,0) self.assertEqual(P._k%1,0)
...@@ -63,14 +63,14 @@ class TestPsfao(unittest.TestCase): ...@@ -63,14 +63,14 @@ class TestPsfao(unittest.TestCase):
def test_bounds_length(self): def test_bounds_length(self):
npix = 100 npix = 100
samp = 0.3 samp = 0.3
P = Psfao((npix,npix),system=load('muse_nfm'),samp=samp) P = Psfao((npix,npix),system=muse_nfm,samp=samp)
self.assertEqual(len(P.bounds[0]),7) self.assertEqual(len(P.bounds[0]),7)
self.assertEqual(len(P.bounds[1]),7) self.assertEqual(len(P.bounds[1]),7)
def test_psd_integral(self): def test_psd_integral(self):
npix = 1024 npix = 1024
samp = 2.0 samp = 2.0
P = Psfao((npix,npix),system=load('zimpol'),samp=samp) P = Psfao((npix,npix),system=zimpol,samp=samp)
fao = P.system.Nact/(2.0*P.system.D) fao = P.system.Nact/(2.0*P.system.D)
df = 1.0/(P.system.D*P._samp_over) df = 1.0/(P.system.D*P._samp_over)
...@@ -104,7 +104,7 @@ class TestPsfao(unittest.TestCase): ...@@ -104,7 +104,7 @@ class TestPsfao(unittest.TestCase):
def test_otf_max(self): def test_otf_max(self):
npix = 1024 npix = 1024
samp = 2.0 samp = 2.0
P = Psfao((npix,npix),system=load('zimpol'),samp=samp) P = Psfao((npix,npix),system=zimpol,samp=samp)
r0 = 0.25 r0 = 0.25
C = 1e-4 C = 1e-4
A = 1.0 A = 1.0
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment