; docformat = 'rst'
;
; NAME: 
;    EVALELLIPSE
; PURPOSE:
;    Evaluate the ordinate (Y) coordinates of points on an ellipse from an array of
;    abscissa (X) and ellpise parameters.
;+
; :Description:
;    Evaluate the ordinate (Y) coordinates of points on an ellipse from an array of
;    abscissa (X) and ellpise parameters.
;
; :Categories:
;    Maths, Evaluation
;    
; :Params:
;    X: in, required, type=float
;        Array (any dimension) of x values.
;    P: in, required, type=float
;        The parameters of the ellipse::
;          P[0] = X center
;          P[1] = Y center
;          P[2] = Semi-major axis
;          P[3] = Semi-minor axis
;          P[4] = Angle between major axis and vertical (deg)
;
; :Returns: The evaluated ordinates at positions X
;
; :Examples:
;    Draw an ellipse of axes 10 and 5, centered on (50,25) and tilted
;    by 20 degrees. Then evaluate the ordinates for integer X.
;      IDL> x=indgen(100)
;      IDL> p=[50,25,10,5,20.]
;      IDL> y= evalellipse(x,p)
;      IDL> ;plot the ellipse
;      IDL> theta=2*!PI*findgen(361)/360.
;      IDL> rx=p(2)*cos(theta)
;      IDL> ry=p(3)*sin(theta)
;      IDL> tilt=(90-p(4))*!DTOR
;      IDL> plot, p(0)+rx*cos(tilt) -ry*sin(tilt), p(1)+rx*sin(tilt) +ry*cos(tilt), /ISO, xr=[0,100]
;      IDL> ;plot the evaluated points 
;      IDL> oplot, x, y(0,*), psym=4
;      IDL> oplot, x, y(1,*), psym=4
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in June 2013, B. Carry (IMCCE)
;      2013 Sep. - B.Carry (IMCCE) - Header for idldoc
;      2017 Jan. - B.Carry (OCA) - idl2 added
;-
function evalEllipse, X, P

 ;--I-- Initialization And Input Verification
 ;--I.1-- Make Compilation Silent
  COMPILE_OPT hidden, idl2

 ;--I.2-- Check Presence of Arguments
  if N_params() LT 2 then begin
    print,'Syntax -  Y =  evalEllipse( X, P )'
    return, -1
 endif

 ;--I.3-- Convert Inputs to Float
  X=float(X)
  P=float(P)

 ;--II-- Variables pre-Computation
  nPoint = n_elements( x )

 ;--II.1-- Center
  x0 = float(p[0])  ;-Center X
  y0 = float(p[1])  ;-Center Y
  xS = x-x0         ;-Shifted X

 ;--II.2-- Radii
  rx = float(p[2])  ;-SMA
  ry = float(p[3])  ;-SmA
  rx2= rx*rx        ;-SMA^2
  ry2= ry*ry        ;-SmA^2

 ;--II.3-- Rotation Parameters
  c = cos(p[4]) ;-Cos
  s = sin(p[4]) ;-Sin
  c2 = c*c      ;-Cos^2
  s2 = s*s      ;-Sin^2
  cs = c*s      ;-Cos*Sin


 ;--III-Discriminant
  dA = replicate( ry2*s2 + rx2*c2, nPoint )
  dB = 2* ( ry2 * ( cs*xS - y0*s2 ) - $
            rx2 * ( cs*xS + y0*c2 )   )
  dC = ry2 * ( xS*xS*c2 + y0*y0*s2 - 2*cs*y0*xS ) + $
       rx2 * ( xS*xS*s2 + y0*y0*c2 + 2*cs*y0*xS ) - $
       rx2 *ry2
  delta = dB*dB - 4*dA*dC

 ;--IV-- Ellipse Evaluation
  return,  yOut = (-[1,1.]#dB + [-1,1.]#sqrt(delta) ) / ([2,2.]#dA)

end