; docformat = 'rst'
;
; NAME: 
;    evalMoffat2d
; PURPOSE:
;    Evaluate the ordinate (Z) coordinates of a 2-D Moffat function
;
;+
; :Description:
;    Evaluate the ordinate (Z) coordinates of a 2-D Moffat function
;
; :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 Moffat 2-D peak::
;          P[0]   Constant baseline level
;          P[1]   Peak value
;          P[2]   Peak half-width (x)
;          P[3]   Peak half-width (y)
;          P[4]   Peak centroid (x)
;          P[5]   Peak centroid (y)
;          P[6]   Rotation angle (radians)
;          P[7]   Moffat power law index - default = 1
;
; :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 August 2017, B. Carry (OCA)
;-
function evalMoffat2d, X, Y, P, _EXTRA=_EXTRA
  COMPILE_OPT hidden, idl2

 ;--------------------------------------------------------------------------------
 ;--I-- Initialization And Input Verification  -----------------------------------
 ;--I.1-- Check Presence of Arguments
  if N_params() LT 3 then begin
    message, /ioError,'Syntax -  Z =  evalMoffat2d( X, Y, P )'
  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-- Moffat Parameters
  widx  = abs(p[2]) > 1e-20 & widy  = abs(p[3]) > 1e-20 
  xp    = xx-p[4]           & yp    = yy-p[5]
  if n_elements(p) ge 8 then power = p[7] else power=1.
  if n_elements(p) ge 7 then theta = p[6] else theta=0.

 ;--II.2-- Moffat Exponent
  if theta NE 0 then begin
    c  = cos(theta) & s  = sin(theta)
    U = ( (xp * (c/widx) - yp * (s/widx))^2 + $
          (xp * (s/widy) + yp * (c/widy))^2 )
  endif else begin
    U = (xp/widx)^2 + (yp/widy)^2
  endelse

 ;--II.2-- 
  mask = u LT (smax^2)  ;; Prevents floating underflow
  moffat = p[1] * mask / (U+1)^power

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

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

 
end