psfmodel.py 25.1 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
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
#%% 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
    """
    
279
    def __init__(self,npix):
FETICK Romain's avatar
FETICK Romain committed
280
281
282
        """
        Parameters
        ----------
283
        npix : tuple of two elements
FETICK Romain's avatar
FETICK Romain committed
284
285
            Model X and Y pixel size when called
        """
286
287
        if type(npix)!=tuple: raise TypeError("Argument `npix` must be a tuple")
        self.npix = npix
rfetick's avatar
rfetick committed
288
        self.bounds = (-_np.inf,_np.inf)
FETICK Romain's avatar
FETICK Romain committed
289
290
    
    def __repr__(self):
291
        return "ParametricPSF of size (%u,%u)"%self.npix
FETICK Romain's avatar
FETICK Romain committed
292
293
294
295
296
297
298
    
    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
299
        return _fft.fftshift(_fft.fft2(psf))
FETICK Romain's avatar
FETICK Romain committed
300
301
302
    
    def mtf(self,*args,**kwargs):
        """Return the Modulation Transfer Function (MTF)"""
rfetick's avatar
rfetick committed
303
        return _np.abs(self.otf(args,kwargs))
FETICK Romain's avatar
FETICK Romain committed
304
305
306
307
    
    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
308
        hdu = _fits.PrimaryHDU(psf, hdr)
FETICK Romain's avatar
FETICK Romain committed
309
310
311
        hdu.writeto(filename, overwrite=True)
        
    def _getfitshdr(self,param,keys=None,keys_comment=None):
312
313
        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`")
FETICK Romain's avatar
FETICK Romain committed
314
        if keys_comment is not None:
315
            if len(keys_comment)!=len(param): raise ValueError("When defined, `keys_comment` must be same size as `param`")
rfetick's avatar
rfetick committed
316
        hdr = _fits.Header()
FETICK Romain's avatar
FETICK Romain committed
317
        
FETICK Romain's avatar
FETICK Romain committed
318
        hdr["HIERARCH ORIGIN"] = "MAOPPY automatic header"
rfetick's avatar
rfetick committed
319
        hdr["HIERARCH CREATION"] = (_time.ctime(),"Date of file creation")
FETICK Romain's avatar
FETICK Romain committed
320
321
322
323
324
325
326
327
328
329
330
331
332
        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
333
        super().__init__(_np.shape(image_psf))
FETICK Romain's avatar
FETICK Romain committed
334
335
336
337
338
339
340
341
        self.image_psf = image_psf
        self.bounds = ()
        
    def __call__(self,*args,**kwargs):
        return self.image_psf


class Moffat(ParametricPSF):
342
343
    def __init__(self,npix,norm=_np.inf):
        self.npix = npix
FETICK Romain's avatar
FETICK Romain committed
344
        self.norm = norm
rfetick's avatar
rfetick committed
345
346
        bounds_down = [_EPSILON,_EPSILON,-_np.inf,1+_EPSILON]
        bounds_up = [_np.inf for i in range(4)]
FETICK Romain's avatar
FETICK Romain committed
347
        self.bounds = (bounds_down,bounds_up)
FETICK Romain's avatar
FETICK Romain committed
348
        
349
350
351
        self._XY = _np.mgrid[0:npix[0],0:npix[1]]
        self._XY[0] -= npix[0]//2
        self._XY[1] -= npix[1]//2
FETICK Romain's avatar
FETICK Romain committed
352
353
354
355
356
357
358
359
360
361
362
    
    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
363
        y = _np.concatenate((x,[dx,dy]))
FETICK Romain's avatar
FETICK Romain committed
364
        return moffat(self._XY,y,norm=self.norm)
FETICK Romain's avatar
FETICK Romain committed
365
366
367
368
369
    
    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
370
class Gaussian(ParametricPSF):
371
372
    def __init__(self,npix):
        self.npix = npix
rfetick's avatar
rfetick committed
373
374
        bounds_down = [_EPSILON,_EPSILON,-_np.inf]
        bounds_up = [_np.inf for i in range(3)]
FETICK Romain's avatar
FETICK Romain committed
375
376
        self.bounds = (bounds_down,bounds_up)
        
377
378
379
        self._XY = _np.mgrid[0:npix[0],0:npix[1]]
        self._XY[1] -= npix[0]/2
        self._XY[0] -= npix[1]/2
FETICK Romain's avatar
FETICK Romain committed
380
381
382
383
384
385
386
387
388
389
    
    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
390
        y = _np.concatenate((x,[dx,dy]))
FETICK Romain's avatar
FETICK Romain committed
391
392
393
        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
394
        super(Gaussian,self).tofits(param,filename,*args,keys=keys,**kwargs)
395

396
397
398
#%% MASTER CLASS
class ParametricPSFfromPSD(ParametricPSF):
    """This class is NOT to be instantiated"""
399
400
401
402
    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")
rfetick's avatar
rfetick committed
403
404
        if system is None: raise ValueError("Keyword `system` must be defined")
        if samp is None: raise ValueError("Keyword `samp` must be defined")
rfetick's avatar
rfetick committed
405
406
        self._npix = npix # "_" to bypass the _xyarray update, that will be made with samp setter
        self.samp = samp # also init _xyarray
rfetick's avatar
rfetick committed
407
        self.system = system
408
409
        self._nparam = nparam
        
rfetick's avatar
rfetick committed
410
411
412
413
414
415
416
417
418
    @property
    def npix(self):
        return self._npix
    
    @npix.setter
    def npix(self,value):
        self._npix = value
        self._computeXYarray()
        
rfetick's avatar
rfetick committed
419
420
421
422
423
424
425
426
    @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)
rfetick's avatar
rfetick committed
427
        self._computeXYarray()
rfetick's avatar
rfetick committed
428
        
rfetick's avatar
rfetick committed
429
430
431
432
433
    def _computeXYarray(self):
        nx,ny = self._nxnyOver
        self._xyarray = _np.mgrid[0:nx, 0:ny].astype(float)
        self._xyarray[0] -= nx//2
        self._xyarray[1] -= ny//2
rfetick's avatar
rfetick committed
434
    
435
    #@_lru_cache(maxsize=2)
436
437
    @property
    def _pix2freq(self):
438
439
440
441
442
        """
        pixel to frequency conversion in the PSD plane
        """
        return 1.0/(self.system.D*self._samp_over)
    
443
    #@_lru_cache(maxsize=2)
444
445
    @property
    def _nxnyOver(self):
446
447
448
        """
        Return the number of pixels for the correctly sampled array (at least at Shannon-Nyquist)
        """
449
        return self.npix[0]*self._k, self.npix[1]*self._k
450
    
451
    #@_lru_cache(maxsize=2)
452
453
454
    @property
    def _nullFreqIndex(self):
        nx,ny = self._nxnyOver
455
456
        return nx//2, ny//2
    
457
    #@_lru_cache(maxsize=2)
458
459
    @property
    def _fxfy(self):
460
461
462
        """
        Return the array of frequencies (correctly sampled above Shannon-Nyquist)
        """
rfetick's avatar
rfetick committed
463
        return self._xyarray * self._pix2freq
464
    
rfetick's avatar
rfetick committed
465
466
467
468
469
470
471
472
473
    def strehlOTF(self,x0):
        otf = self.otf(x0,_caller='self')
        otf_diff = self._otf_diffraction()
        return _np.real(_np.sum(otf)/_np.sum(otf_diff))
    
    def strehlMarechal(self,x0):
        _,sig2 = self.psd(x0)
        return _np.exp(-sig2)
    
474
475
476
    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
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    def check_parameters(self,x0):
        bd, bu = self.bounds
        x0 = _np.array(x0)
        bd = _np.array(bd)
        bu = _np.array(bu)
        if len(x0)!=self._nparam: raise ValueError('len(x0) is different from length of bounds')
        if _np.any(x0<bd): raise ValueError('Lower bounds are not respected')
        if _np.any(x0>bu): raise ValueError('Upper bounds are not respected')
        
    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'):
            raise ValueError("Cannot call `otf(...)` when undersampled (functionality not implemented yet)")
        
        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
    
rfetick's avatar
rfetick committed
503
    def _otf_turbulent(self,x0):
504
        PSD,integral = self.psd(x0)
rfetick's avatar
rfetick committed
505
        L = self.system.D * self._samp_over
rfetick's avatar
rfetick committed
506
        Bg = _fft.fft2(_fft.fftshift(PSD)) / L**2
507
508
509
        #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
510
    
511
    #@_lru_cache(maxsize=2)
rfetick's avatar
rfetick committed
512
    def _otf_diffraction(self):
rfetick's avatar
rfetick committed
513
514
515
516
        nx, ny = self._nxnyOver
        NpupX = _np.ceil(nx/self._samp_over)
        NpupY = _np.ceil(ny/self._samp_over)
        tab = _np.zeros((nx, ny), dtype=_np.complex)
rfetick's avatar
rfetick committed
517
        tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
rfetick's avatar
rfetick committed
518
        return _fft.fftshift(_np.abs(_fft.ifft2(_np.abs(_fft.fft2(tab)) ** 2)) / _np.sum(tab))
rfetick's avatar
rfetick committed
519
    
520
    #@_lru_cache(maxsize=2)
rfetick's avatar
rfetick committed
521
    def _otf_shift(self,dx,dy):
rfetick's avatar
rfetick committed
522
523
524
        Y, X = self._xyarray
        ny,nx = X.shape
        return _np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny))
rfetick's avatar
rfetick committed
525
526
527
528
529
530
    
    def __call__(self,x0,dx=0,dy=0):
        """
        Parameters
        ----------
        x0 : numpy.array (dim=1), tuple, list
531
            Array of parameters for the PSD (see __doc__)
rfetick's avatar
rfetick committed
532
533
534
535
536
537
        dx : float
            PSF X shifting [pix] (default = 0)
        dy : float
            PSF Y shifting [pix] (default = 0)
        """
        
rfetick's avatar
rfetick committed
538
        out = _np.real(_fft.fftshift(_fft.ifft2(_fft.fftshift(self.otf(x0,dx=dx,dy=dy,_caller='self')))))
rfetick's avatar
rfetick committed
539
540
541
542
        
        if self._k==1:
            return out
        else:
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
            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]
    """
559
    def __init__(self,npix,system=None,samp=None):
560
561
562
        """
        Parameters
        ----------
563
        npix : tuple
564
565
566
567
568
569
            Size of output PSF
        system : OpticalSystem
            Optical system for this PSF
        samp : float
            Sampling at the observation wavelength
        """
570
        super().__init__(2,npix,system=system,samp=samp)
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
        
        # 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)
593
594
        nx0,ny0 = self._nullFreqIndex
        f2D = self._fxfy
595
596
597
598
599
600
601
602
603
604
        
        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
605
606
        fmax = _np.min([nx0,ny0])*self._pix2freq
        integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum
607
608
609
610
        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
        integral = integral_in + integral_out
        
        return PSD, integral
611

612
#%% PSFAO MODEL        
613
class Psfao(ParametricPSFfromPSD):
FETICK Romain's avatar
FETICK Romain committed
614
    """
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
    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
634
    """
635
    def __init__(self,npix,system=None,Lext=10.,samp=None):
FETICK Romain's avatar
FETICK Romain committed
636
637
638
        """
        Parameters
        ----------
639
        npix : tuple
FETICK Romain's avatar
FETICK Romain committed
640
641
642
643
644
645
646
647
648
            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
        """
649
        super().__init__(7,npix,system=system,samp=samp)
650
        
FETICK Romain's avatar
FETICK Romain committed
651
        self.Lext = Lext
rfetick's avatar
rfetick committed
652
        # r0,C,A,alpha,ratio,theta,beta 
rfetick's avatar
rfetick committed
653
654
        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
655
656
657
658
659
660
661
662
663
        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
664
            See __doc__ for more details
FETICK Romain's avatar
FETICK Romain committed
665
666
667
668
            
        Returns
        -------
        psd : numpy.array (dim=2)
669
        integral : float : the integral of the `psd` up to infinity
FETICK Romain's avatar
FETICK Romain committed
670
671
            
        """
672
        self.check_parameters(x0)
673
674
        nx0,ny0 = self._nullFreqIndex
        f2D = self._fxfy
FETICK Romain's avatar
FETICK Romain committed
675
        
676
677
        F2 = f2D[0]**2. + f2D[1]**2.
        
FETICK Romain's avatar
FETICK Romain committed
678
679
        Fao = self.system.Nact/(2.0*self.system.D)
        
rfetick's avatar
rfetick committed
680
        r0,C,A,alpha,ratio,theta,beta = x0
681
682
        
        PSD = 0.0229* r0**(-5./3.) * ((1. / self.Lext**2.) + F2)**(-11./6.)
FETICK Romain's avatar
FETICK Romain committed
683
684
        PSD *= (F2 >= Fao**2.)
        
rfetick's avatar
rfetick committed
685
686
        ax = alpha*ratio
        ay = alpha/ratio
687
        param = (ax,ay,theta,beta,0,0)
688
689
690
691
        
        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
rfetick's avatar
rfetick committed
692
        # Warning: Moffat numerical normalization generates strehlOTF jump when "_k" is changed
693
694
695
        
        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
696
697
        
        # Compute PSD integral up to infinity
698
699
        fmax = _np.min([nx0,ny0])*self._pix2freq
        integral_in = _np.sum(PSD*(F2<(fmax**2))) * self._pix2freq**2 # numerical sum
rfetick's avatar
rfetick committed
700
        integral_out = 0.0229*6*_np.pi/5 * (r0*fmax)**(-5./3.) # analytical sum
701
702
703
        integral = integral_in + integral_out
        
        return PSD, integral
FETICK Romain's avatar
FETICK Romain committed
704
        
rfetick's avatar
rfetick committed
705
    def tofits(self,param,filename,*args,keys=None,overwrite=False,**kwargs):
FETICK Romain's avatar
FETICK Romain committed
706
        if keys is None:
rfetick's avatar
rfetick committed
707
            keys = ["R0","CST","SIGMA2","ALPHA","RATIO","THETA","BETA"]
FETICK Romain's avatar
FETICK Romain committed
708
709
710
            keys_comment = ["Fried parameter [m]",
                            "PSD AO area constant C [rad2]",
                            "PSD AO area Moffat variance A [rad2]",
711
                            "PSD AO area Moffat alpha [1/m]",
rfetick's avatar
rfetick committed
712
                            "PSD AO area Moffat sqrt(ax/ay) ratio",
FETICK Romain's avatar
FETICK Romain committed
713
714
                            "PSD AO area Moffat theta [rad]",
                            "PSD AO area Moffat beta"]
FETICK Romain's avatar
FETICK Romain committed
715
716
717
718
719
720
        
        # 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")
721
722
        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
723
724
725
        hdr["HIERARCH SAMP"] = (self.samp,"Sampling (eg. 2 for Shannon)")
        hdr["HIERARCH LEXT"] = (self.Lext,"Von-Karman outer scale")
        
rfetick's avatar
rfetick committed
726
        hdu = _fits.PrimaryHDU(psf, hdr)
rfetick's avatar
rfetick committed
727
        hdu.writeto(filename, overwrite=overwrite)