Commit 750942b7 authored by skamann's avatar skamann
Browse files

Modified how FWHM is estimated for MAOPPY profile.

parent a47aa9e1
......@@ -457,7 +457,7 @@ class Parameters(object):
def polynomial_fit(self, name, order, mask=None, polynom="plain"):
def polynomial_fit(self, name, order, mask=None, polynom="plain", bounds=None):
Fit a parameter using a polynomial.
......@@ -478,6 +478,13 @@ class Parameters(object):
polynom : str, optional
The function used in the fitting. May be either 'plain',
'legendre', 'hermite', or 'chebyshev'.
bounds : tuple
In case the polynomial fit should be
success : bool
Flag indicating if the fit has been successful.
# sanity check on provided parameter
if name not in self.names:
......@@ -503,8 +510,11 @@ class Parameters(object):
mask = self.mask[(name, 'value')]
# carry out fit
fit =[~mask],[~mask, (name, 'value')], order)[(name, 'fit')] = fit(abscissa)
f =[~mask],[~mask, (name, 'value')], order)
fitted_values = f(abscissa)[(name, 'fit')] = fitted_values
def make_hdu(self, include_fit=True, header_keys=None):
......@@ -35,12 +35,13 @@ import io
import logging
import numpy as np
from contextlib import redirect_stdout
from scipy.interpolate import RectBivariateSpline, interp2d
from scipy.interpolate import RectBivariateSpline, InterpolatedUnivariateSpline, interp2d
from .profile import Profile, FitResult
if importlib.util.find_spec('maoppy') is not None:
from maoppy.psfmodel import Psfao, psffit, oversample
from maoppy.instrument import muse_nfm
from maoppy.utils import circavgplt
logging.error('Package "maoppy" is not installed. Cannot use MAOPPY PSF model.')
......@@ -424,7 +425,6 @@ class Maoppy(Profile):
centre : tuple of 2 floats
The central coordinates of the PSF along the y- and x-axes.
if shape is None:
shape = self.maoppy_shape
return shape[0]/2 - (k - 1.)/(2.*k), shape[1]/2 - (k - 1.)/(2.*k)
......@@ -432,9 +432,14 @@ class Maoppy(Profile):
def calculate_fwhm(cls, r0, bck, amp, alpha, e, beta, wvl=6e-7, **kwargs):
Calculate and return the FWHM of the MAOPPY PSF profile for a given
Estimate and return the FWHM of the MAOPPY PSF profile for a given
set of parameters.
To estimate the FWHM, the method creates the MAOPPY PSF on a 2dim.
grid, then determines the symmetric radial profile, and tries to
obtain the radii where the intensity drops to half its central value
using univariate spline interpolation.
r0 : float
......@@ -469,11 +474,18 @@ class Maoppy(Profile):
raise IOError('Unknown argument(s) to call of Maoppy.fwhm: {0}'.format(kwargs))
psf = Psfao(Npix=(30, 30), system=muse_nfm, samp=muse_nfm.samp(wvl))
theta = 0.
f = psf((r0, bck, amp, alpha, 1-e, theta, beta))[:15]
f = psf((r0, bck, amp, alpha, 1 - e, theta, beta))
r, fr = circavgplt(f)
f = InterpolatedUnivariateSpline(r, fr - f.max() / 2.)
roots = f.roots()
return 2.*np.sqrt(np.sum(f > (0.5*f.max()))/np.pi)
if len(roots) != 2:
logger.error("Could not estimate FWHM of MAOPPY profile.")
return np.nan
return roots.max() - roots.min()
def fwhm(self):
......@@ -483,16 +495,8 @@ class Maoppy(Profile):
fwhm : float
The FWHM of the MAOPPY PSF profile with the given parameters.
from scipy.stats import binned_statistic
# get mean radii for fluxes for annuli of widths ~1 pixel
bstat = binned_statistic(self.r, [self.r, self.f], bins=self.maxrad)
r_bin, f_bin = bstat.statistic
# interpolate along flux axis to get FWHM. Note that flux values must be monotonically INCREASING
fwhm = 2.*np.interp(0.5*self.f.max(), f_bin[::-1], r_bin[::-1])
return fwhm
return Maoppy.calculate_fwhm(r0=self.r0, bck=self.bck, amp=self.amp, alpha=self.alpha, e=self.e,
beta=self.beta, wvl=self.wvl)
def strehl(self):
......@@ -1433,6 +1433,14 @@ class Profile(Canvas):
return radii, np.cumsum(cog_flux)
def calculate_fwhm(cls, **kwargs):
Calculate FWHM for a given set of parameters. Must be overwritten by
classes inheriting from Profile().
raise NotImplementedError
def fwhm(self):
......@@ -263,7 +263,8 @@ class Variables(Parameters):
self.lut.index = index
if not self.lut.index.isin(index).all():
logger.warning('Changing dispersion axis as requested will alter existing values of look-up table.')
'Changing dispersion axis as requested will alter existing values of look-up table.')
self.lut = self.lut.reindex(index=index, method='nearest')
def fit_lut(self, order, mask=None):
Markdown is supported
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