; docformat = 'rst'
;
; NAME: 
;    evalPSFAO
; PURPOSE:
;    Evaluate the Point-Spread Function of an Adaptative-opics fed Camera
;
;+
; :Description:
;    Evaluate the Point-Spread Function of an Adaptative-opics fed Camera
;
; :Categories:
;    Maths, Evaluation
;    
; :Params:
;    X: in, required, type=float
;        Array (any dimension) of x values.
;    Y: in, required, type=float
;        Array (any dimension) of y values.
;    P: in, required, type=float
;        The parameters of the Gaussian 2d peak::
;          P[0]   Constant baseline level
;          P[1]   Core Peak value
;          P[2]   Core Peak half-width (x)
;          P[3]   Core Peak half-width (y)
;          P[4]   Core Peak centroid (x)
;          P[5]   Core Peak centroid (y)
;          P[6]   Core Rotation angle (radians)
;          P[7]   Halo Peak value
;          P[8]   Halo Peak half-width (x)
;          P[9]   Halo Peak half-width (y)
;          P[10]  Halo Peak centroid (x)
;          P[11]  Halo Peak centroid (y)
;          P[12]  Halo Rotation angle (radians)
;
; :Returns: An array of [Nx,Ny] elements with evaluated ordinates at positions X,Y.
;
; :Keywords:
;    _EXTRA: in, optional, type=struct
;        Can suppleant P in inversion code s(see MPFITFUN)
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in April 2015, B. Carry (IMCCE)
;      2017 Aug. - B.Carry (OCA) - Added idl2 compile option - Core is now a Lorentzian function
;-
function evalPSFAO, X, Y, P, _EXTRA=_EXTRA
  COMPILE_OPT hidden, idl2
 ;--------------------------------------------------------------------------------
 ;--I-- Initialization And Input Verification  -----------------------------------
 ;--I.1-- Check Presence of Arguments
  if N_params() LT 2 then begin
    print,'Syntax -  Z =  evalPSFAO( X, Y, P )'
    return, -1
  endif

 ;--I.2-- Handle Extra Keyword (e.g., for use in MPFIT): replace P by _EXTRA.P
  if keyword_set( _extra ) then p = _extra.p

 ;--I.3-- Convert 1-D Inputs to 2-D Arrays
  dimX=size(X)
  dimY=size(Y)
  if dimX[0] NE 2 then xx = x[*] # (y[*]*0 + 1)   else xx = x
  if dimY[0] NE 2 then yy = (x[*]*0 + 1) # y[*]   else yy = y

  P=float(P)
  sz = size(xx)
  if sz[sz[0]+1] EQ 5 then smax = 26D else smax = 13.


 ;--------------------------------------------------------------------------------
 ;--II-- Function Evaluation  ----------------------------------------------------
 ;--II.1-- Core PSF
  core = evalMoffat2d( x, y, [0,p[1:6],1] )

 ;--II.2-- Halo Seeing
  halo = evalGauss2d( x, y, [p[0],p[7:12]] )


 ;--------------------------------------------------------------------------------
 ;--III-- Background Evaluation  -------------------------------------------------
  plane= p[0]


 ;--------------------------------------------------------------------------------
 ;--IV-- Return Gaussian 2D Peak with Background  --------------------------------
  return, (core+halo+plane)

 
end