psfmodel.py 24.3 KB
Newer Older
FETICK Romain's avatar
FETICK Romain committed
1
2
3
4
5
6
7
8
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon May 27 17:31:18 2019

@author: rfetick
"""

rfetick's avatar
rfetick committed
9
10
11
12
import numpy as _np
from scipy.optimize import least_squares as _least_squares
from astropy.io import fits as _fits
import time as _time
rfetick's avatar
rfetick committed
13
import numpy.fft as _fft # import fft2, fftshift, ifft2
rfetick's avatar
rfetick committed
14
15
import sys as _sys
from maoppy.utils import binning as _binning
rfetick's avatar
rfetick committed
16
from functools import lru_cache as _lru_cache
FETICK Romain's avatar
FETICK Romain committed
17

18
# Value to compute finite differences
rfetick's avatar
rfetick committed
19
_EPSILON = _np.sqrt(_sys.float_info.epsilon)
20

rfetick's avatar
rfetick committed
21
22
23
#%% FUNCTIONS
def oversample(samp):
    """Find the minimal integer that allows oversampling"""
rfetick's avatar
rfetick committed
24
    k = int(_np.ceil(2.0/samp))
rfetick's avatar
rfetick committed
25
26
27
    return (k*samp,k)

#%% FITTING
FETICK Romain's avatar
FETICK Romain committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def lsq_flux_bck(model, data, weights, background=True, positive_bck=False):
    """Compute the analytical least-square solution for flux and background
    LS = SUM_pixels { weights*(flux*model + bck - data)² }
    
    Parameters
    ----------
    model: numpy.ndarray
    data: numpy.ndarray
    weights: numpy.ndarray
    
    Keywords
    --------
    background: bool
        Activate/inactivate background (activated by default:True)
    positive_bck : bool
        Makes background positive (default:False)
    """
rfetick's avatar
rfetick committed
45
46
47
48
49
    ws = _np.sum(weights)
    mws = _np.sum(model * weights)
    mwds = _np.sum(model * weights * data)
    m2ws = _np.sum(weights * (model ** 2))
    wds = _np.sum(weights * data)
FETICK Romain's avatar
FETICK Romain committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

    if background:
        delta = mws ** 2 - ws * m2ws
        amp = 1. / delta * (mws * wds - ws * mwds)
        bck = 1. / delta * (-m2ws * wds + mws * mwds)
    else:
        amp = mwds / m2ws
        bck = 0.0
        
    if bck<0 and positive_bck: #re-implement above equation
        amp = mwds / m2ws
        bck = 0.0

    return amp, bck

def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
66
           positive_bck=False,fixed=None,npixfit=None,**kwargs):
FETICK Romain's avatar
FETICK Romain committed
67
68
69
70
71
72
73
74
75
76
77
78
79
    """Fit a PSF with a parametric model solving the least-square problem
       epsilon(x) = SUM_pixel { weights * (amp * Model(x) + bck - psf)² }
    
    Parameters
    ----------
    psf : numpy.ndarray
        The experimental image to be fitted
    Model : class
        The class representing the fitting model
    x0 : tuple, list, numpy.ndarray
        Initial guess for parameters
    weights : numpy.ndarray
        Least-square weighting matrix (same size as `psf`)
FETICK Romain's avatar
FETICK Romain committed
80
        Inverse of the noise's variance
FETICK Romain's avatar
FETICK Romain committed
81
82
83
84
85
86
87
88
89
90
        Default: uniform weighting
    dxdy : tuple of two floats
        Eventual guess on PSF shifting
    flux_bck : tuple of two bool
        Only background can be activate/inactivated
        Flux is always activated (sorry!)
    positive_bck : bool
        Force background to be positive or null
    fixed : numpy.ndarray
        Fix some parameters to their initial value (default: None)
91
92
    npixfit : int (default=None)
        Increased pixel size for improved fitting accuracy
FETICK Romain's avatar
FETICK Romain committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
    **kwargs :
        All keywords used to instantiate your `Model`
    
    Returns
    -------
    out.x : numpy.array
            Parameters at optimum
       .dxdy : tuple of 2 floats
           PSF shift at optimum
       .flux_bck : tuple of two floats
           Estimated image flux and background
       .psf : numpy.ndarray (dim=2)
           Image of the PSF model at optimum
       .success : bool
           Minimization success
       .status : int
           Minimization status (see scipy doc)
       .message : string
           Human readable minimization status
       .active_mask : numpy.array
           Saturated bounds
       .nfev : int
           Number of function evaluations
       .cost : float
           Value of cost function at optimum
    """
119
    
rfetick's avatar
rfetick committed
120
    if weights is None: weights = _np.ones_like(psf)
121
    elif len(psf)!=len(weights): raise ValueError("Keyword `weights` must have same number of elements as `psf`")
122
123
124
    
    # Increase array size for improved fitting accuracy
    if npixfit is not None:
rfetick's avatar
rfetick committed
125
        sx,sy = _np.shape(psf)
126
        if (sx>npixfit) or (sy>npixfit): raise ValueError('npixfit must be greater or equal to both psf dimensions')
rfetick's avatar
rfetick committed
127
128
        psfbig = _np.zeros((npixfit,npixfit))
        wbig = _np.zeros((npixfit,npixfit))
129
130
131
132
133
        psfbig[npixfit//2-sx//2:npixfit//2+sx//2,npixfit//2-sy//2:npixfit//2+sy//2] = psf
        wbig[npixfit//2-sx//2:npixfit//2+sx//2,npixfit//2-sy//2:npixfit//2+sy//2] = weights
        psf = psfbig
        weights = wbig
    
rfetick's avatar
rfetick committed
134
135
    sqW = _np.sqrt(weights)
    Model_inst = Model(_np.shape(psf),**kwargs)
FETICK Romain's avatar
FETICK Romain committed
136
137
138
139
140
    
    class CostClass(object):
        def __init__(self):
            self.iter = 0
        def __call__(self,y):
141
            if (self.iter%3)==0: print("-",end="")
FETICK Romain's avatar
FETICK Romain committed
142
143
144
145
            self.iter += 1
            x, dxdy = mini2input(y)
            m = Model_inst(x,dx=dxdy[0],dy=dxdy[1])
            amp, bck = lsq_flux_bck(m, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
rfetick's avatar
rfetick committed
146
            return _np.reshape(sqW * (amp * m + bck - psf), _np.size(psf))
FETICK Romain's avatar
FETICK Romain committed
147
148
149
150
    
    cost = CostClass()
    
    if fixed is not None:
151
        if len(fixed)!=len(x0): raise ValueError("When defined, `fixed` must be same size as `x0`")
FETICK Romain's avatar
FETICK Romain committed
152
        FREE = [not fixed[i] for i in range(len(fixed))]
rfetick's avatar
rfetick committed
153
        INDFREE = _np.where(FREE)[0]
FETICK Romain's avatar
FETICK Romain committed
154
155
156
    
    def input2mini(x,dxdy):
        # Transform user parameters to minimizer parameters
157
        if fixed is None: xfree = x
rfetick's avatar
rfetick committed
158
159
        else: xfree = _np.take(x,INDFREE)
        return _np.concatenate((xfree,dxdy))
FETICK Romain's avatar
FETICK Romain committed
160
161
162
163
164
165
    
    def mini2input(y):
        # Transform minimizer parameters to user parameters
        if fixed is None:
            xall = y[0:-2]
        else:
rfetick's avatar
rfetick committed
166
            xall = _np.copy(x0)
FETICK Romain's avatar
FETICK Romain committed
167
168
169
170
171
172
            for i in range(len(y)-2):
                xall[INDFREE[i]] = y[i]
        return (xall,y[-2:])
    
    def get_bound(inst):
        b_low = inst.bounds[0]
rfetick's avatar
rfetick committed
173
174
        if fixed is not None: b_low = _np.take(b_low,INDFREE)
        b_low = _np.concatenate((b_low,[-_np.inf,-_np.inf]))
FETICK Romain's avatar
FETICK Romain committed
175
        b_up = inst.bounds[1]
rfetick's avatar
rfetick committed
176
177
        if fixed is not None: b_up = _np.take(b_up,INDFREE)
        b_up = _np.concatenate((b_up,[_np.inf,_np.inf]))
FETICK Romain's avatar
FETICK Romain committed
178
179
        return (b_low,b_up)
    
rfetick's avatar
rfetick committed
180
    result = _least_squares(cost, input2mini(x0,dxdy), bounds=get_bound(Model_inst))
FETICK Romain's avatar
FETICK Romain committed
181
182
183
184
185
186
187
188
189
190
191
192
    
    print("") #finish line of "-"
    
    result.x, result.dxdy = mini2input(result.x) #split output between x and dxdy

    m = Model_inst(result.x,dx=result.dxdy[0],dy=result.dxdy[1])
    amp, bck = lsq_flux_bck(m, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
    
    result.flux_bck = (amp,bck)
    result.psf = m    
    return result

193
#%% BASIC FUNCTIONS
194
def reduced_coord(XY,ax,ay,theta,cx,cy):
rfetick's avatar
rfetick committed
195
196
197
    c = _np.cos(theta)
    s = _np.sin(theta)
    s2 = _np.sin(2.0 * theta)
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225

    Rxx = (c/ax)**2 + (s/ay)**2
    Ryy = (c/ay)**2 + (s/ax)**2
    Rxy =  s2/ay**2 -  s2/ax**2
    
    u = Rxx*(XY[0]-cx)**2 + Rxy*(XY[0]-cx)*(XY[1]-cy) + Ryy*(XY[1]-cy)**2
    return u    

def moffat(XY, param, norm=None):
    """
    Compute a Moffat function on a meshgrid
    moff = Enorm * (1+u)^(-beta)
    with `u` the quadratic coordinates in the shifted and rotated frame
    
    Parameters
    ----------
    XY : numpy.ndarray (dim=2)
        The (X,Y) meshgrid with X = XY[0] and Y = XY[1]
    param : list, tuple, numpy.ndarray (len=6)
        param[0] - Alpha X
        param[1] - Alpha Y
        param[2] - Theta
        param[3] - Beta
        param[4] - Center X
        param[5] - Center Y
        
    Keywords
    --------
rfetick's avatar
rfetick committed
226
    norm : None, _np.inf, float (>0), int (>0)
227
228
229
        Radius for energy normalization
        None      - No energy normalization (maximum=1.0)
                    Enorm = 1.0
rfetick's avatar
rfetick committed
230
        _np.inf    - Total energy normalization (on the whole X-Y plane)
231
232
233
234
235
236
237
238
239
240
241
                    Enorm = (beta-1)/(pi*ax*ay)
        float,int - Energy normalization up to the radius defined by this value
                    Enorm = (beta-1)/(pi*ax*ay)*(1-(1+(R**2)/(ax*ay))**(1-beta))    
    """
    if len(param)!=6: raise ValueError("Parameter `param` must contain exactly 6 elements, but input has %u elements"%len(param))
    ax,ay,theta,beta,cx,cy = param
    
    u = reduced_coord(XY,ax,ay,theta,cx,cy)
    
    if norm is None:
        Enorm = 1
rfetick's avatar
rfetick committed
242
    elif norm == _np.inf:
243
        if beta<=1: raise ValueError("Cannot compute Moffat energy for beta<=1")
rfetick's avatar
rfetick committed
244
        Enorm = (beta-1) / (_np.pi*ax*ay)
245
246
    else:
        if beta==1: raise ValueError("Energy computation for beta=1.0 not implemented yet. Sorry!")
rfetick's avatar
rfetick committed
247
        Enorm = (beta-1) / (_np.pi*ax*ay)
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
        Enorm = Enorm / (1 - (1 + (norm**2) / (ax*ay))**(1-beta))
    
    return Enorm * (1.0+u)**(-beta)

def gauss(XY,param):
    """
    Compute a Gaussian function on a meshgrid

    Parameters
    ----------
    XY : numpy.ndarray (dim=2)
        The (X,Y) meshgrid with X = XY[0] and Y = XY[1]
    param : list, tuple, numpy.ndarray (len=5)
        param[0] - Sigma X
        param[1] - Sigma Y
        param[2] - Theta
        param[3] - Center X
        param[4] - Center Y  
    """
    if len(param)!=5: raise ValueError("Parameter `param` must contain exactly 5 elements, but input has %u elements"%len(param))
rfetick's avatar
rfetick committed
268
269
    ax = _np.sqrt(2)*param[0]
    ay = _np.sqrt(2)*param[1]
270
    u = reduced_coord(XY,ax,ay,param[2],param[3],param[4])
rfetick's avatar
rfetick committed
271
    return 1.0/(2*_np.pi*param[0]*param[1])*_np.exp(-u)
272

FETICK Romain's avatar
FETICK Romain committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
#%% CLASS PARAMETRIC PSF AND ITS SUBCLASSES
class ParametricPSF(object):
    """Super-class defining parametric PSFs
    Not to be instantiated, only serves as a referent for subclasses
    """
    
    def __init__(self,Npix):
        """
        Parameters
        ----------
        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
rfetick's avatar
rfetick committed
289
        self.bounds = (-_np.inf,_np.inf)
FETICK Romain's avatar
FETICK Romain committed
290
291
292
293
294
295
296
297
298
299
    
    def __repr__(self):
        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")
    
    def otf(self,*args,**kwargs):
        """Return the Optical Transfer Function (OTF)"""
        psf = self.__call__(args,kwargs)
rfetick's avatar
rfetick committed
300
        return _fft.fftshift(_fft.fft2(psf))
FETICK Romain's avatar
FETICK Romain committed
301
302
303
    
    def mtf(self,*args,**kwargs):
        """Return the Modulation Transfer Function (MTF)"""
rfetick's avatar
rfetick committed
304
        return _np.abs(self.otf(args,kwargs))
FETICK Romain's avatar
FETICK Romain committed
305
306
307
308
    
    def tofits(self,param,filename,*args,keys=None,keys_comment=None,**kwargs):
        psf = self.__call__(param,*args,**kwargs)
        hdr = self._getfitshdr(param,keys=keys,keys_comment=keys_comment)
rfetick's avatar
rfetick committed
309
        hdu = _fits.PrimaryHDU(psf, hdr)
FETICK Romain's avatar
FETICK Romain committed
310
311
312
313
314
315
316
317
318
319
        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_comment is not None:
            if len(keys_comment)!=len(param):
                raise ValueError("When defined, `keys_comment` must be same size as `param`")
rfetick's avatar
rfetick committed
320
        hdr = _fits.Header()
FETICK Romain's avatar
FETICK Romain committed
321
        
FETICK Romain's avatar
FETICK Romain committed
322
        hdr["HIERARCH ORIGIN"] = "MAOPPY automatic header"
rfetick's avatar
rfetick committed
323
        hdr["HIERARCH CREATION"] = (_time.ctime(),"Date of file creation")
FETICK Romain's avatar
FETICK Romain committed
324
325
326
327
328
329
330
331
332
333
334
335
336
        for i in range(len(param)):
            if keys_comment is None:
                hdr["HIERARCH PARAM "+keys[i]] = param[i]
            else:
                hdr["HIERARCH PARAM "+keys[i]] = (param[i],keys_comment[i])
        return hdr


class ConstantPSF(ParametricPSF):
    """Create a constant PSF, given as a 2D image, using ParametricPSF formalism
    With such a formalism, a constant PSF is just a particular case of a parametric PSF
    """
    def __init__(self,image_psf):
rfetick's avatar
rfetick committed
337
        super().__init__(_np.shape(image_psf))
FETICK Romain's avatar
FETICK Romain committed
338
339
340
341
342
343
344
345
        self.image_psf = image_psf
        self.bounds = ()
        
    def __call__(self,*args,**kwargs):
        return self.image_psf


class Moffat(ParametricPSF):
rfetick's avatar
rfetick committed
346
    def __init__(self,Npix,norm=_np.inf):
FETICK Romain's avatar
FETICK Romain committed
347
348
        self.Npix = Npix
        self.norm = norm
rfetick's avatar
rfetick committed
349
350
        bounds_down = [_EPSILON,_EPSILON,-_np.inf,1+_EPSILON]
        bounds_up = [_np.inf for i in range(4)]
FETICK Romain's avatar
FETICK Romain committed
351
        self.bounds = (bounds_down,bounds_up)
FETICK Romain's avatar
FETICK Romain committed
352
        
rfetick's avatar
rfetick committed
353
        self._XY = _np.mgrid[0:Npix[0],0:Npix[1]]
FETICK Romain's avatar
FETICK Romain committed
354
355
        self._XY[1] -= Npix[0]/2
        self._XY[0] -= Npix[1]/2
FETICK Romain's avatar
FETICK Romain committed
356
357
358
359
360
361
362
363
364
365
366
    
    def __call__(self,x,dx=0,dy=0):
        """
        Parameters
        ----------
        x : list, tuple, numpy.ndarray (len=4)
            x[0] - Alpha X
            x[1] - Alpha Y
            x[2] - Theta
            x[3] - Beta
        """
rfetick's avatar
rfetick committed
367
        y = _np.concatenate((x,[dx,dy]))
FETICK Romain's avatar
FETICK Romain committed
368
        return moffat(self._XY,y,norm=self.norm)
FETICK Romain's avatar
FETICK Romain committed
369
370
371
372
373
    
    def tofits(self,param,filename,*args,keys=["ALPHA_X","ALPHA_Y","THETA","BETA"],**kwargs):
        super(Moffat,self).tofits(param,filename,*args,keys=keys,**kwargs)


FETICK Romain's avatar
FETICK Romain committed
374
375
376
class Gaussian(ParametricPSF):
    def __init__(self,Npix):
        self.Npix = Npix
rfetick's avatar
rfetick committed
377
378
        bounds_down = [_EPSILON,_EPSILON,-_np.inf]
        bounds_up = [_np.inf for i in range(3)]
FETICK Romain's avatar
FETICK Romain committed
379
380
        self.bounds = (bounds_down,bounds_up)
        
rfetick's avatar
rfetick committed
381
        self._XY = _np.mgrid[0:Npix[0],0:Npix[1]]
FETICK Romain's avatar
FETICK Romain committed
382
383
384
385
386
387
388
389
390
391
392
393
        self._XY[1] -= Npix[0]/2
        self._XY[0] -= Npix[1]/2
    
    def __call__(self,x,dx=0,dy=0):
        """
        Parameters
        ----------
        x : list, tuple, numpy.ndarray (len=4)
            x[0] - Sigma X
            x[1] - Sigma Y
            x[2] - Theta
        """
rfetick's avatar
rfetick committed
394
        y = _np.concatenate((x,[dx,dy]))
FETICK Romain's avatar
FETICK Romain committed
395
396
397
        return gauss(self._XY,y)
    
    def tofits(self,param,filename,*args,keys=["SIGMA_X","SIGMA_Y","THETA"],**kwargs):
FETICK Romain's avatar
FETICK Romain committed
398
        super(Gaussian,self).tofits(param,filename,*args,keys=keys,**kwargs)
399

400
401
402
403
#%% MASTER CLASS
class ParametricPSFfromPSD(ParametricPSF):
    """This class is NOT to be instantiated"""
    def __init__(self,nparam,Npix,system=None,samp=None):
rfetick's avatar
rfetick committed
404
        if not (type(Npix) in [tuple,list,_np.ndarray]): raise ValueError("Npix must be a tuple, list or numpy.ndarray")
rfetick's avatar
rfetick committed
405
406
407
408
409
410
411
        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.system = system
        self.samp = samp
412
413
        self._nparam = nparam
        
rfetick's avatar
rfetick committed
414
415
416
417
418
419
420
421
    @property
    def samp(self):
        return self._samp_over/self._k
    
    @samp.setter
    def samp(self,value):
        # Manage cases of undersampling
        self._samp_over, self._k = oversample(value)
422

rfetick's avatar
rfetick committed
423
424
    def check_parameters(self,x0):
        bd, bu = self.bounds
rfetick's avatar
rfetick committed
425
426
427
        x0 = _np.array(x0)
        bd = _np.array(bd)
        bu = _np.array(bu)
428
        if len(x0)!=self._nparam: raise ValueError('len(x0) is different from length of bounds')
rfetick's avatar
rfetick committed
429
430
        if _np.any(x0<bd): raise ValueError('Lower bounds are not respected')
        if _np.any(x0>bu): raise ValueError('Upper bounds are not respected')
rfetick's avatar
rfetick committed
431
432
433
434
435
436
437
438
439
440
        
    def otf(self,x0,dx=0,dy=0,_caller='user'):
        """
        See __call__ for input arguments
        Warning: result of otf will be unconsistent if undersampled!!!
        This issue is solved with oversampling + binning in __call__ but not here
        For the moment, the `_caller` keyword prevents user to misuse otf
        """
        
        if (self._k > 1) and (_caller != 'self'):
441
            raise ValueError("Cannot call `otf(...)` when undersampled (functionality not implemented yet)")
rfetick's avatar
rfetick committed
442
443
444
445
446
447
448
        
        OTF_TURBULENT = self._otf_turbulent(x0)
        OTF_DIFFRACTION = self._otf_diffraction()
        OTF_SHIFT = self._otf_shift(dx,dy)
        
        return OTF_TURBULENT * OTF_DIFFRACTION * OTF_SHIFT
    
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
    def getpix2freq(self):
        """
        pixel to frequency conversion in the PSD plane
        """
        return 1.0/(self.system.D*self._samp_over)
    
    def getnxny_over(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
    
    def get_null_freq_index(self):
        nx,ny = self.getnxny_over()
        return nx//2, ny//2
    
    def getfxfy(self):
        """
        Return the array of frequencies (correctly sampled above Shannon-Nyquist)
        """
        Nx_over, Ny_over = self.getnxny_over()
        f2D = _np.mgrid[0:Nx_over, 0:Ny_over].astype(float)
        f2D[0] -= Nx_over//2
        f2D[1] -= Ny_over//2
        f2D *= self.getpix2freq()
        return f2D
    
    def psd(self,x0):
        raise ValueError("ParametricPSFfromPSD is not to be instantiated. the `psd` method must be override in the subclasses")
    
rfetick's avatar
rfetick committed
479
    def _otf_turbulent(self,x0):
480
        PSD,integral = self.psd(x0)
rfetick's avatar
rfetick committed
481
        L = self.system.D * self._samp_over
rfetick's avatar
rfetick committed
482
        Bg = _fft.fft2(_fft.fftshift(PSD)) / L**2
483
484
485
        #Dphi = _fft.fftshift(_np.real(2 * (Bg[0, 0] - Bg))) # normalized on the numerical FoV 
        Dphi = _fft.fftshift(_np.real(2 * (integral - Bg))) # normalized up to infinity
        return _np.exp(-Dphi/2.)
rfetick's avatar
rfetick committed
486
    
rfetick's avatar
rfetick committed
487
    @_lru_cache(maxsize=2)
rfetick's avatar
rfetick committed
488
489
490
491
    def _otf_diffraction(self):
        Nx_over = self.Npix[0]*self._k
        Ny_over = self.Npix[1]*self._k
        
rfetick's avatar
rfetick committed
492
493
494
        NpupX = _np.ceil(Nx_over/self._samp_over)
        NpupY = _np.ceil(Ny_over/self._samp_over)
        tab = _np.zeros((Nx_over, Ny_over), dtype=_np.complex)
rfetick's avatar
rfetick committed
495
        tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
rfetick's avatar
rfetick committed
496
        return _fft.fftshift(_np.abs(_fft.ifft2(_np.abs(_fft.fft2(tab)) ** 2)) / _np.sum(tab))
rfetick's avatar
rfetick committed
497
498
499
500
501
    
    def _otf_shift(self,dx,dy):
        Nx_over = self.Npix[0]*self._k
        Ny_over = self.Npix[1]*self._k
        
rfetick's avatar
rfetick committed
502
503
504
505
        Y, X = _np.mgrid[0:Nx_over,0:Ny_over].astype(float)
        X = (X-Nx_over/2) * 2*_np.pi*1j/Nx_over
        Y = (Y-Ny_over/2) * 2*_np.pi*1j/Ny_over
        return _np.exp(-X*dx*self._k - Y*dy*self._k)
rfetick's avatar
rfetick committed
506
507
508
509
510
511
    
    def __call__(self,x0,dx=0,dy=0):
        """
        Parameters
        ----------
        x0 : numpy.array (dim=1), tuple, list
512
            Array of parameters for the PSD (see __doc__)
rfetick's avatar
rfetick committed
513
514
515
516
517
518
        dx : float
            PSF X shifting [pix] (default = 0)
        dy : float
            PSF Y shifting [pix] (default = 0)
        """
        
rfetick's avatar
rfetick committed
519
        out = _np.real(_fft.fftshift(_fft.ifft2(_fft.fftshift(self.otf(x0,dx=dx,dy=dy,_caller='self')))))
rfetick's avatar
rfetick committed
520
521
522
523
        
        if self._k==1:
            return out
        else:
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
            return _binning(out,int(self._k)) # undersample PSF if needed (if it was oversampled for computations)

#%% TURBULENT PSF
class Turbulent(ParametricPSFfromPSD):
    """
    Summary
    -------
    PSF model dedicated to long-exposure imaging with turbulence
    p = Turbulent((npix,npix),system=system,samp=samp)
    
    Description
    -----------
    the PSF is parametrised through the PSD of the electromagnetic phase
        x[0] - Fried parameter r0 [m]
        x[1] - Von Karman external length [m]
    """
    def __init__(self,Npix,system=None,samp=None):
        """
        Parameters
        ----------
        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)
        
        # r0,L0
        bounds_down = [_EPSILON, _EPSILON]
        bounds_up = [_np.inf, _np.inf]
        self.bounds = (bounds_down,bounds_up)
    
    def psd(self,x0):
        """Compute the PSD model from parameters
        PSD is given in [rad²/f²] = [rad² m²]
        
        Parameters
        ----------
        x0 : numpy.array (dim=1), tuple, list
            See __doc__ for more details
            
        Returns
        -------
        psd : numpy.array (dim=2)
        integral : float : the integral of the `psd` up to infinity
            
        """
        self.check_parameters(x0)
        nx0,ny0 = self.get_null_freq_index()
        f2D = self.getfxfy()
        pix2freq = self.getpix2freq()
        
        F2 = f2D[0]**2. + f2D[1]**2.
        
        r0,Lext = x0
        PSD = 0.0229* r0**(-5./3.) * ((1./Lext**2.) + F2)**(-11./6.)

        # Set PSD = 0 at null frequency
        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
        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
        integral = integral_in + integral_out
        
        return PSD, integral
593

rfetick's avatar
rfetick committed
594
#%% PSFAO MODEL
595
class Psfao(ParametricPSFfromPSD):
FETICK Romain's avatar
FETICK Romain committed
596
    """
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
    Summary
    -------
    PSF model dedicated to long-exposure imaging with adaptive optics
    p = Psfao((npix,npix),system=system,samp=samp)
    
    Description
    -----------
    the PSF is parametrised through the PSD of the electromagnetic phase
        x[0] - Fried parameter           : r0 [m]
        x[1] - Corrected area background : C [rad² m²]
        x[2] - Moffat phase variance     : A [rad²]
        x[3] - Moffat width              : alpha [1/m]
        x[4] - Moffat elongation ratio   : sqrt(ax/ay)
        x[5] - Moffat angle              : theta [rad]
        x[6] - Moffat power law          : beta
        
    Reference
    ---------
    Fétick et al., August 2019, A&A, Vol.628
FETICK Romain's avatar
FETICK Romain committed
616
    """
617
    def __init__(self,Npix,system=None,Lext=10.,samp=None):
FETICK Romain's avatar
FETICK Romain committed
618
619
620
621
622
623
624
625
626
627
628
629
630
        """
        Parameters
        ----------
        Npix : tuple
            Size of output PSF
        system : OpticalSystem
            Optical system for this PSF
        samp : float
            Sampling at the observation wavelength
        Lext : float
            Von-Karman external scale (default = 10 m)
            Useless if Fao >> 1/Lext
        """
631
        super().__init__(7,Npix,system=system,samp=samp)
632
        
FETICK Romain's avatar
FETICK Romain committed
633
        self.Lext = Lext
rfetick's avatar
rfetick committed
634
        # r0,C,A,alpha,ratio,theta,beta 
rfetick's avatar
rfetick committed
635
636
        bounds_down = [_EPSILON,0,0,_EPSILON,_EPSILON,-_np.inf,1+_EPSILON]
        bounds_up = [_np.inf for i in range(7)]
FETICK Romain's avatar
FETICK Romain committed
637
638
639
640
641
642
643
644
645
        self.bounds = (bounds_down,bounds_up)
    
    def psd(self,x0):
        """Compute the PSD model from parameters
        PSD is given in [rad²/f²] = [rad² m²]
        
        Parameters
        ----------
        x0 : numpy.array (dim=1), tuple, list
646
            See __doc__ for more details
FETICK Romain's avatar
FETICK Romain committed
647
648
649
650
            
        Returns
        -------
        psd : numpy.array (dim=2)
651
        integral : float : the integral of the `psd` up to infinity
FETICK Romain's avatar
FETICK Romain committed
652
653
            
        """
654
        self.check_parameters(x0)
655
656
657
        nx0,ny0 = self.get_null_freq_index()
        f2D = self.getfxfy()
        pix2freq = self.getpix2freq()
FETICK Romain's avatar
FETICK Romain committed
658
        
659
660
        F2 = f2D[0]**2. + f2D[1]**2.
        
FETICK Romain's avatar
FETICK Romain committed
661
662
        Fao = self.system.Nact/(2.0*self.system.D)
        
rfetick's avatar
rfetick committed
663
        r0,C,A,alpha,ratio,theta,beta = x0
664
665
        
        PSD = 0.0229* r0**(-5./3.) * ((1. / self.Lext**2.) + F2)**(-11./6.)
FETICK Romain's avatar
FETICK Romain committed
666
667
        PSD *= (F2 >= Fao**2.)
        
rfetick's avatar
rfetick committed
668
669
        ax = alpha*ratio
        ay = alpha/ratio
670
        param = (ax,ay,theta,beta,0,0)
rfetick's avatar
rfetick committed
671
        PSD += (F2 < Fao**2.) * _np.abs(C + A*moffat(f2D,param,norm=Fao))
672
        # Set PSD = 0 at null frequency
673
        PSD[nx0,ny0] = 0.0
674
675
        
        # Compute PSD integral up to infinity
676
        fmax = _np.min([nx0,ny0])*pix2freq
rfetick's avatar
rfetick committed
677
678
        integral_in = _np.sum(PSD*(F2<(fmax**2))) * pix2freq**2 # numerical sum
        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
679
680
681
        integral = integral_in + integral_out
        
        return PSD, integral
FETICK Romain's avatar
FETICK Romain committed
682
        
rfetick's avatar
rfetick committed
683
    def tofits(self,param,filename,*args,keys=None,overwrite=False,**kwargs):
FETICK Romain's avatar
FETICK Romain committed
684
        if keys is None:
rfetick's avatar
rfetick committed
685
            keys = ["R0","CST","SIGMA2","ALPHA","RATIO","THETA","BETA"]
FETICK Romain's avatar
FETICK Romain committed
686
687
688
            keys_comment = ["Fried parameter [m]",
                            "PSD AO area constant C [rad2]",
                            "PSD AO area Moffat variance A [rad2]",
689
                            "PSD AO area Moffat alpha [1/m]",
rfetick's avatar
rfetick committed
690
                            "PSD AO area Moffat sqrt(ax/ay) ratio",
FETICK Romain's avatar
FETICK Romain committed
691
692
                            "PSD AO area Moffat theta [rad]",
                            "PSD AO area Moffat beta"]
FETICK Romain's avatar
FETICK Romain committed
693
694
695
696
697
698
        
        # redefine tofits() because extra hdr is required
        psf = self.__call__(param,*args,**kwargs)
        hdr = self._getfitshdr(param,keys=keys,keys_comment=keys_comment)
        
        hdr["HIERARCH SYSTEM"] = (self.system.name,"System name")
699
700
        hdr["HIERARCH SYSTEM D"] = (self.system.D,"Primary mirror diameter")
        hdr["HIERARCH SYSTEM NACT"] = (self.system.Nact,"Linear number of AO actuators")
FETICK Romain's avatar
FETICK Romain committed
701
702
703
        hdr["HIERARCH SAMP"] = (self.samp,"Sampling (eg. 2 for Shannon)")
        hdr["HIERARCH LEXT"] = (self.Lext,"Von-Karman outer scale")
        
rfetick's avatar
rfetick committed
704
        hdu = _fits.PrimaryHDU(psf, hdr)
rfetick's avatar
rfetick committed
705
        hdu.writeto(filename, overwrite=overwrite)