Commit f7a8f1ad authored by rfetick's avatar rfetick
Browse files

apply oversampling shift compensation internally

parent a0834cfb
......@@ -46,8 +46,8 @@ class PsfaoGalpak(Psfao, PointSpreadFunction):
ww = wvl_um/self.wvl_ref_um # wavelength ratio
r0,C,A,alpha,ratio,theta,beta = self.param
p = [r0*ww**(6./5.),C*ww**(-2),A*ww**(-2),alpha,ratio,theta,beta]
dx = (self._k-1)/(2*self._k) # centered on one pixel
dy = (self._k-1)/(2*self._k) # centered on one pixel
dx = 0
dy = 0
return self.__call__(p, dx=dx, dy=dy)
......
......@@ -10,10 +10,9 @@ import numpy as _np
from scipy.optimize import least_squares as _least_squares
from astropy.io import fits as _fits
import time as _time
import numpy.fft as _fft # import fft2, fftshift, ifft2
import numpy.fft as _fft
import sys as _sys
from maoppy.utils import binning as _binning
#from functools import lru_cache as _lru_cache
# Value to compute finite differences
_EPSILON = _np.sqrt(_sys.float_info.epsilon)
......@@ -74,6 +73,7 @@ def lsq_flux_bck(model, data, weights, background=True, positive_bck=False):
return amp, bck
def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
positive_bck=False,fixed=None,npixfit=None,**kwargs):
"""Fit a PSF with a parametric model solving the least-square problem
......@@ -202,6 +202,7 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
result.psf = m
return result
#%% BASIC FUNCTIONS
def reduced_coord(XY,ax,ay,theta,cx,cy,grad=False):
"""Create an array of reduced coordinated with elongation and rotation"""
......@@ -234,6 +235,7 @@ def reduced_coord(XY,ax,ay,theta,cx,cy,grad=False):
return u
def moffat(XY, param, norm=None, removeInside=0, grad=False):
"""Compute a Moffat function on a meshgrid
```moff = E * (1+u)^(-beta)```
......@@ -307,6 +309,7 @@ def moffat(XY, param, norm=None, removeInside=0, grad=False):
return V*E*F
def gauss(XY,param):
"""
Compute a Gaussian function on a meshgrid
......@@ -328,6 +331,7 @@ def gauss(XY,param):
u = reduced_coord(XY,ax,ay,param[2],param[3],param[4])
return 1.0/(2*_np.pi*param[0]*param[1])*_np.exp(-u)
#%% CLASS PARAMETRIC PSF AND ITS SUBCLASSES
class ParametricPSF(object):
"""Super-class defining parametric PSFs
......@@ -451,6 +455,7 @@ class Gaussian(ParametricPSF):
def tofits(self,param,filename,*args,keys=["SIGMA_X","SIGMA_Y","THETA"],**kwargs):
super(Gaussian,self).tofits(param,filename,*args,keys=keys,**kwargs)
#%% MASTER CLASS
class ParametricPSFfromPSD(ParametricPSF):
"""This class is NOT to be instantiated"""
......@@ -493,7 +498,6 @@ class ParametricPSFfromPSD(ParametricPSF):
self._xyarray[0] -= nx//2
self._xyarray[1] -= ny//2
#@_lru_cache(maxsize=2)
@property
def _pix2freq(self):
"""
......@@ -501,7 +505,6 @@ class ParametricPSFfromPSD(ParametricPSF):
"""
return 1.0/(self.system.D*self._samp_over)
#@_lru_cache(maxsize=2)
@property
def _nxnyOver(self):
"""
......@@ -509,13 +512,11 @@ class ParametricPSFfromPSD(ParametricPSF):
"""
return self.npix[0]*self._k, self.npix[1]*self._k
#@_lru_cache(maxsize=2)
@property
def _nullFreqIndex(self):
nx,ny = self._nxnyOver
return nx//2, ny//2
#@_lru_cache(maxsize=2)
@property
def _fxfy(self):
"""
......@@ -591,7 +592,6 @@ class ParametricPSFfromPSD(ParametricPSF):
return otf, g2
return otf
#@_lru_cache(maxsize=2)
def _otfDiffraction(self):
nx, ny = self._nxnyOver
NpupX = _np.ceil(nx/self._samp_over)
......@@ -600,10 +600,12 @@ class ParametricPSFfromPSD(ParametricPSF):
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
return _fft.fftshift(_np.abs(_fft.ifft2(_np.abs(_fft.fft2(tab)) ** 2)) / _np.sum(tab))
#@_lru_cache(maxsize=2)
def _otfShift(self,dx,dy):
Y, X = self._xyarray
ny,nx = X.shape
# Compensate oversampling shift
dx += (self._k-1)/(2*self._k)
dy += (self._k-1)/(2*self._k)
return _np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny))
def __call__(self,x0,dx=0,dy=0,grad=False):
......@@ -639,6 +641,7 @@ class ParametricPSFfromPSD(ParametricPSF):
return _binning(psf,k), g2
return _binning(psf,k) # undersample PSF if needed (if it was oversampled for computations)
#%% TURBULENT PSF
class Turbulent(ParametricPSFfromPSD):
"""
......@@ -717,6 +720,7 @@ class Turbulent(ParametricPSFfromPSD):
if grad: return PSD,integral,g,integral_g
return PSD, integral
#%% PSFAO MODEL
class Psfao(ParametricPSFfromPSD):
"""
......
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