Commit e032360c authored by rfetick's avatar rfetick
Browse files

Psfao moffat is numerically normalized, modified some internal syntaxes

parent f46feb9e
......@@ -276,20 +276,19 @@ class ParametricPSF(object):
Not to be instantiated, only serves as a referent for subclasses
"""
def __init__(self,Npix):
def __init__(self,npix):
"""
Parameters
----------
Npix : tuple of two elements
npix : tuple of two elements
Model X and Y pixel size when called
"""
if type(Npix)!=tuple:
raise TypeError("Argument `Npix` must be a tuple")
self.Npix = Npix
if type(npix)!=tuple: raise TypeError("Argument `npix` must be a tuple")
self.npix = npix
self.bounds = (-_np.inf,_np.inf)
def __repr__(self):
return "ParametricPSF of size (%u,%u)"%self.Npix
return "ParametricPSF of size (%u,%u)"%self.npix
def __call__(self,*args,**kwargs):
raise ValueError("ParametricPSF is not made to be instantiated. Better use its subclasses")
......@@ -310,13 +309,10 @@ class ParametricPSF(object):
hdu.writeto(filename, overwrite=True)
def _getfitshdr(self,param,keys=None,keys_comment=None):
if keys is None:
keys = ["PARAM %u"%(i+1) for i in range(len(param))]
if len(keys)!=len(param):
raise ValueError("When defined, `keys` must be same size as `param`")
if keys is None: keys = ["PARAM %u"%(i+1) for i in range(len(param))]
if len(keys)!=len(param): raise ValueError("When defined, `keys` must be same size as `param`")
if keys_comment is not None:
if len(keys_comment)!=len(param):
raise ValueError("When defined, `keys_comment` must be same size as `param`")
if len(keys_comment)!=len(param): raise ValueError("When defined, `keys_comment` must be same size as `param`")
hdr = _fits.Header()
hdr["HIERARCH ORIGIN"] = "MAOPPY automatic header"
......@@ -343,16 +339,16 @@ class ConstantPSF(ParametricPSF):
class Moffat(ParametricPSF):
def __init__(self,Npix,norm=_np.inf):
self.Npix = Npix
def __init__(self,npix,norm=_np.inf):
self.npix = npix
self.norm = norm
bounds_down = [_EPSILON,_EPSILON,-_np.inf,1+_EPSILON]
bounds_up = [_np.inf for i in range(4)]
self.bounds = (bounds_down,bounds_up)
self._XY = _np.mgrid[0:Npix[0],0:Npix[1]]
self._XY[0] -= Npix[0]//2
self._XY[1] -= Npix[1]//2
self._XY = _np.mgrid[0:npix[0],0:npix[1]]
self._XY[0] -= npix[0]//2
self._XY[1] -= npix[1]//2
def __call__(self,x,dx=0,dy=0):
"""
......@@ -372,15 +368,15 @@ class Moffat(ParametricPSF):
class Gaussian(ParametricPSF):
def __init__(self,Npix):
self.Npix = Npix
def __init__(self,npix):
self.npix = npix
bounds_down = [_EPSILON,_EPSILON,-_np.inf]
bounds_up = [_np.inf for i in range(3)]
self.bounds = (bounds_down,bounds_up)
self._XY = _np.mgrid[0:Npix[0],0:Npix[1]]
self._XY[1] -= Npix[0]/2
self._XY[0] -= Npix[1]/2
self._XY = _np.mgrid[0:npix[0],0:npix[1]]
self._XY[1] -= npix[0]/2
self._XY[0] -= npix[1]/2
def __call__(self,x,dx=0,dy=0):
"""
......@@ -400,13 +396,13 @@ class Gaussian(ParametricPSF):
#%% MASTER CLASS
class ParametricPSFfromPSD(ParametricPSF):
"""This class is NOT to be instantiated"""
def __init__(self,nparam,Npix,system=None,samp=None):
if not (type(Npix) in [tuple,list,_np.ndarray]): raise ValueError("Npix must be a tuple, list or numpy.ndarray")
if len(Npix)!=2: raise ValueError("Npix must be of length = 2")
if (Npix[0]%2) or (Npix[1]%2): raise ValueError("Each Npix component must be even")
def __init__(self,nparam,npix,system=None,samp=None):
if not (type(npix) in [tuple,list,_np.ndarray]): raise ValueError("npix must be a tuple, list or numpy.ndarray")
if len(npix)!=2: raise ValueError("npix must be of length = 2")
if (npix[0]%2) or (npix[1]%2): raise ValueError("Each npix component must be even")
if system is None: raise ValueError("Keyword `system` must be defined")
if samp is None: raise ValueError("Keyword `samp` must be defined")
self.Npix = Npix
self.npix = npix
self.system = system
self.samp = samp
self._nparam = nparam
......@@ -447,34 +443,38 @@ class ParametricPSFfromPSD(ParametricPSF):
return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
#@_lru_cache(maxsize=2)
def getpix2freq(self):
@property
def _pix2freq(self):
"""
pixel to frequency conversion in the PSD plane
"""
return 1.0/(self.system.D*self._samp_over)
#@_lru_cache(maxsize=2)
def getnxny_over(self):
@property
def _nxnyOver(self):
"""
Return the number of pixels for the correctly sampled array (at least at Shannon-Nyquist)
"""
return self.Npix[0]*self._k, self.Npix[1]*self._k
return self.npix[0]*self._k, self.npix[1]*self._k
#@_lru_cache(maxsize=2)
def get_null_freq_index(self):
nx,ny = self.getnxny_over()
@property
def _nullFreqIndex(self):
nx,ny = self._nxnyOver
return nx//2, ny//2
#@_lru_cache(maxsize=2)
def getfxfy(self):
@property
def _fxfy(self):
"""
Return the array of frequencies (correctly sampled above Shannon-Nyquist)
"""
Nx_over, Ny_over = self.getnxny_over()
Nx_over, Ny_over = self._nxnyOver
f2D = _np.mgrid[0:Nx_over, 0:Ny_over].astype(float)
f2D[0] -= Nx_over//2
f2D[1] -= Ny_over//2
f2D *= self.getpix2freq()
f2D *= self._pix2freq
return f2D
def strehlOTF(self,x0):
......@@ -499,8 +499,8 @@ class ParametricPSFfromPSD(ParametricPSF):
#@_lru_cache(maxsize=2)
def _otf_diffraction(self):
Nx_over = self.Npix[0]*self._k
Ny_over = self.Npix[1]*self._k
Nx_over = self.npix[0]*self._k
Ny_over = self.npix[1]*self._k
NpupX = _np.ceil(Nx_over/self._samp_over)
NpupY = _np.ceil(Ny_over/self._samp_over)
......@@ -510,8 +510,8 @@ class ParametricPSFfromPSD(ParametricPSF):
#@_lru_cache(maxsize=2)
def _otf_shift(self,dx,dy):
Nx_over = self.Npix[0]*self._k
Ny_over = self.Npix[1]*self._k
Nx_over = self.npix[0]*self._k
Ny_over = self.npix[1]*self._k
Y, X = _np.mgrid[0:Nx_over,0:Ny_over].astype(float)
X = (X-Nx_over/2) * 2*_np.pi*1j/Nx_over
......@@ -551,18 +551,18 @@ class Turbulent(ParametricPSFfromPSD):
x[0] - Fried parameter r0 [m]
x[1] - Von Karman external length [m]
"""
def __init__(self,Npix,system=None,samp=None):
def __init__(self,npix,system=None,samp=None):
"""
Parameters
----------
Npix : tuple
npix : tuple
Size of output PSF
system : OpticalSystem
Optical system for this PSF
samp : float
Sampling at the observation wavelength
"""
super().__init__(2,Npix,system=system,samp=samp)
super().__init__(2,npix,system=system,samp=samp)
# r0,L0
bounds_down = [_EPSILON, _EPSILON]
......@@ -585,9 +585,8 @@ class Turbulent(ParametricPSFfromPSD):
"""
self.check_parameters(x0)
nx0,ny0 = self.get_null_freq_index()
f2D = self.getfxfy()
pix2freq = self.getpix2freq()
nx0,ny0 = self._nullFreqIndex
f2D = self._fxfy
F2 = f2D[0]**2. + f2D[1]**2.
......@@ -598,8 +597,8 @@ class Turbulent(ParametricPSFfromPSD):
PSD[nx0,ny0] = 0.0
# Compute PSD integral up to infinity
fmax = _np.min([nx0,ny0])*pix2freq
integral_in = _np.sum(PSD*(F2<(fmax**2))) * pix2freq**2 # numerical sum
fmax = _np.min([nx0,ny0])*self._pix2freq
integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum
integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
integral = integral_in + integral_out
......@@ -628,11 +627,11 @@ class Psfao(ParametricPSFfromPSD):
---------
Fétick et al., August 2019, A&A, Vol.628
"""
def __init__(self,Npix,system=None,Lext=10.,samp=None):
def __init__(self,npix,system=None,Lext=10.,samp=None):
"""
Parameters
----------
Npix : tuple
npix : tuple
Size of output PSF
system : OpticalSystem
Optical system for this PSF
......@@ -642,7 +641,7 @@ class Psfao(ParametricPSFfromPSD):
Von-Karman external scale (default = 10 m)
Useless if Fao >> 1/Lext
"""
super().__init__(7,Npix,system=system,samp=samp)
super().__init__(7,npix,system=system,samp=samp)
self.Lext = Lext
# r0,C,A,alpha,ratio,theta,beta
......@@ -666,9 +665,8 @@ class Psfao(ParametricPSFfromPSD):
"""
self.check_parameters(x0)
nx0,ny0 = self.get_null_freq_index()
f2D = self.getfxfy()
pix2freq = self.getpix2freq()
nx0,ny0 = self._nullFreqIndex
f2D = self._fxfy
F2 = f2D[0]**2. + f2D[1]**2.
......@@ -682,13 +680,17 @@ class Psfao(ParametricPSFfromPSD):
ax = alpha*ratio
ay = alpha/ratio
param = (ax,ay,theta,beta,0,0)
PSD += (F2 < Fao**2.) * _np.abs(C + A*moffat(f2D,param,norm=Fao))
# Set PSD = 0 at null frequency
PSD[nx0,ny0] = 0.0
moff = moffat(f2D,param,norm=Fao) * (F2 < Fao**2.)
moff[nx0,ny0] = 0.0 # Set Moffat PSD = 0 at null frequency
moff = moff / (_np.sum(moff)*self._pix2freq**2) # normalize moffat numerically to get correct A=sigma² in the AO zone
PSD += (F2 < Fao**2.) * _np.abs(C + A*moff) # before I wrote "moffat(f2D,param,norm=Fao)" instead of "moff", so no numerical normalization
PSD[nx0,ny0] = 0.0 # Set PSD = 0 at null frequency
# Compute PSD integral up to infinity
fmax = _np.min([nx0,ny0])*pix2freq
integral_in = _np.sum(PSD*(F2<(fmax**2))) * pix2freq**2 # numerical sum
fmax = _np.min([nx0,ny0])*self._pix2freq
integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum
integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
integral = integral_in + integral_out
......
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